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PLANAR  ARRAYS  ON  LATTICES  AND  THEIR  FFT  STEERING,  A  PRIMER 


1.  INTRODUCTION 

A  2007  NRL  memo  report  |[Q  formulated  the  problem  of  steering  a  narrowband  planar  receive  array  to 
many  regularly  spaced  beam  directions  simultaneously  using  a  generalized  two-dimensional  (2D)  discrete 
Fourier  transform  (DFT)  and  then  detailed  the  design  of  custom  fast  Fourier  transform  (FFT)  algorithms 
for  realizing  such  steering  digitally.  The  mathematical  heart  of  the  development  was  the  notion  of  modular 
arithmetic  in  2D  point  lattices,  of  taking  a  lattice  modulo  one  of  its  sublattices  or,  equivalently,  of  taking 
an  integer  vector  modulo  a  nonsingular  square  integer  matrix.  These  are  standard  notions  from  elementary 
group  theory,  so  the  report  was  largely  tutorial  and  simply  aimed  to  give  the  engineer  with  no  background 
in  abstract  algebra  a  way  to  acquire  the  tools  needed  to  work  with  these  algorithms. 

The  subsequent  development  of  example  designs  of  such  array-steering  systems  for  a  new,  second 
document — it  will  be  released  eventually  as  a  formal  NRL  report  0 — led  to  an  approach  to  the  design  of 
the  necessary  custom  FFTs  that  reduced  the  fancy  modular  mathematics  to  the  geometry  of  various  tilings 
of  the  plane.  While  this  approach  is  far  more  engineer  friendly,  the  inclusion  of  a  fully  adequate  tutorial  on 
the  tiling-based  approach  would  have  made  that  document  overlong  and  entirely  unwieldy. 

So  a  third  document,  this  report,  was  called  for.  The  primary  topic  here  is  the  tiling  approach.  However, 
in  the  spirit  of  enabling  one-stop  background  shopping  prior  to  tackling  either  of  the  first  two  documents  or 
our  group’s  other  research  papers  on  arrays,  substantial  additional  material  is  included:  (1)  a  basic  array- 
output  formulation  using  minimal  electromagnetic  theory,  (2)  those  ideas  from  multidimensional  digital 
signal  processing  (DSP)  that  are  important  in  array-pattern  and  array-steering  work,  (3)  such  basics  of  the 
theory  of  point  lattices  as  are  needed  to  connect  the  DSP  and  array  ideas,  and  (4)  even  some  key  underlying 
ideas  from  linear  algebra  that  are  perhaps  less  well  grasped  by  many  engineers  than  would  be  ideal.  There 
are  but  two  design  examples  in  this  document,  and  they  aim  not  so  much  to  demonstrate  the  range  of  what 
the  FFT-steering  approach  can  handle  as  to  demonstrate  that  the  simple  cases  can  be  handled  simply.  More 
extensive  design  examples  will  appear  in  the  second  document. 

While  it  is  certainly  reasonable  then  to  view  this  third  report  as  a  gentle,  more  manageable  introduction  to 
the  second,  it  is  just  as  reasonable  to  view  it  as  a  general-purpose  tutorial,  aimed  at  the  nonspecialist  electrical 
or  even  computer  engineer  with  no  prior  array  background,  on  basic  array  ideas  and  on  the  notations  used 
to  represent  them  in  a  whole  family  of  array  papers.  Little  is  assumed  beyond  some  comfort  with  matrix 
algebra  and  familiarity  with  basic  signals  and  systems,  particularly  those  aspects  concerned  with  elementary 
DSP.  The  one  exception  is  the  Appendix,  which  shows  how  the  generalized  2D  DFT  in  certain  special  cases 
can  be  realized  as  an  ordinary  ID  DFT  simply  through  correct  bin  labeling.  That  Appendix  was  left  a  touch 
more  mathematically  demanding  simply  to  keep  it  relatively  brief. 


Manuscript  approved  March  4,  201 1 . 


1 


2 


Jeffrey  O.  Coleman 


1.1  Structure  of  this  Report  and  its  Relation  to  the  Literature 

The  usual  time  and  frequency  domains  of  signals  and  systems  are  dual  in  the  sense  that 


•  functions  of  one  domain  are  Fourier  transformed,  or  inverse  transformed,  to  functions  of  the  other, 
and 

•  when  time  t  and  frequency  f  are  swapped  in  a  Fourier-transform  property,  we  speak  of  the  new  prop¬ 
erty  obtained  (which  might  require  a  sign  change  or  a  factor  of  27t  somewhere)  as  the  dual  of  the 
first. 


Of  course  we  might  use  a  discrete-time  index  n  to  represent  time,  but  the  essence  of  the  matter  is  the  same. 

Likewise,  planar-array  work  involves  a  position  domain  and  a  beamspace  domain  that  arc  dual  in  very 
much  the  same  way.  Each  of  them  is  two  dimensional,  technically  in  the  plane  of  the  array  in  fact,  and 
the  first  order  of  business  below,  in  Section  [2j  is  to  review  enough  linear  algebra  and  point-lattice  basics  to 
enable  us  to  represent  points  of  interest  in  those  two  domains  in  a  convenient  way.  Using  the  lattice  material 
just  developed,  Section[3]then  presents  a  signal-processing  engineer’s  simplified,  minimal-electromagnetics 
view  of  the  simplest  of  planar  arrays:  those  that  are  periodic.  Section |4]re turns  to  tool  building  and  presents 
the  basics  of  working  with  DSP  in  the  position  <->  beamspace  Fourier  domains,  including  the  basics  of  a 
2D  discrete-space  Fourier  transform,  which  turns  out  to  be  nothing  at  all  mysterious  but  simply  an  ordinary 
discrete-time  Fourier  transform  used  twice,  once  for  each  spatial  dimension.  In  Section |5}  the  heart  of  this 
document,  the  material  of  the  preceding  sections  finally  comes  together  in  an  array  theory  simple  enough 
to  be  written  in  a  few  lines  but  detailed  enough  to  permit  us  to  design,  as  the  heart  of  that  section,  a  simple 
FFT  steering  system  for  an  all-digital  periodic  planar  array. 

The  approach  to  the  2D  FFT  of  Section  [5]  is  in  fact  the  unified  Cooley-Tukey  approach  presented  by 
Mersereau  and  Speake  H  in  1981,  though  this  development  differs  markedly  from  theirs  and  aims  to  be  both 
simpler  and  much  more  application  oriented.  This  report  should,  however,  adequately  prepare  the  reader  to 
tackle  that  classic  paper  and  so  gain  a  sense  of  how  these  algorithms  are  dealt  with  in  the  multidimensional- 
DSP  community.  In  applications  where  the  twiddle-factor  multiplications  of  the  Cooley-Tukey  approach 
look  to  be  an  undue  computational  burden,  the  alternative  FFT  approach  of  Guessoum  and  Mersereau  Q  is 
certainly  worth  a  look,  as  it  uses  the  Smith  normal  form  of  an  integer  matrix  to  avoid  the  need  for  twiddle 
factors.  The  absence  of  those  multiplications  helps  computational  efficiency,  but  the  approach  is  far  less 
straightforward  to  grasp. 


Every  FFT  algorithm  discussed  above  implements  a  generalized  2D  DFT  in  which  the  familiar  number 
of  points  N  of  an  ordinary  ID  DFT  has  been  replaced  with  some  2x2  nonsingular  integer  matrix  R  so  that 
1/N  in  the  exponent  of  the  complex  exponential  becomes  R  ,  and  so  forth.  Factorization  of  this  integer 
size  parameter  |[0  or  “periodicity  matrix’  ’  mm  r  into  smaller  matrices  of  the  same  type,  where  “small” 
refers  to  |R|,  the  absolute  value  of  the  determinant,  is  what  leads  to  the  FFT.  As  a  matter  of  algebra,  when 
|R|  is  prime,  such  factorization  is  not  possible  and  the  2D  FFT  reduces  to  a  straightforward,  brute-force  2D 
DFT  calculation.  The  surprise  is  that,  even  though  the  generalized  DFT  expression  shows  a  factor  of  R  1  in 
the  exponent  of  the  complex  exponential,  for  prime  |R|  that  generalized  DFT  is  just  an  ordinary  ID  DFT  of 


|R|  points  with  carefully  chosen  2D  bin  labels.  This  is  proved  in  an  Appendix 
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1.2  Notation 

The  standard  mathematical  notation  for  the  set  of  integers,  all  of  them  including  the  negatives  and  zero, 
is  Z  and  comes  from  the  German  verb  zalilen,  to  count.  It  is  customary  to  write  the  set  of  all  possible  N- 
vectors  with  elements  from  Z  as  ZiV,  but  some  authors  use  this  notation  to  mean  column  vectors,  and  others 
use  it  to  mean  rows.  This  document  uses  both  integer  row  vectors  and  integer  column  vectors  a  great  deal, 
so  the  usual  7LN notation  is  modified  here:  the  set  of  integer  column  two- vectors  is  '%?,  and  the  set  of  integer 
row  two-vectors  is  Z2. 

Other  than  the  above  and  in  complex  exponentials  e.i27r(whatcvcb  powers  arise  rarely  in  this  document 
and  are  applied  directly  only  to  quantities  in  mathematical  “containers”  like  (  ),  |  |,  or  ||  ||.  This  lets  us  use 
superscripts  relatively  freely  for  identifiers.  We  can  speak  of  signal  s1,  signal  s2,  etc.,  with  the  latter  being 
“s  two”  and  not  “s  squared.”  Subscripts  are  frequently  used  as  identifiers  also,  but  most  often  with  some 
spatial  meaning.  Just  as  subscripts  identify  the  position  of  an  element  in  a  vector  or  matrix  so  that  n  =  [  "I,  ] 
or  R  =  [  rr'2\  rll  ]  >  here  subscripts  can  indicate  position  within  some  signal  that  has  elements  spread  over 
space  or  time,  so  that  s2  is  the  spatial  sample  of  signal  s2  at  the  position  indicated  by  n. 

Table  [T]  summarizes  this  document’s  notational  conventions  for  quantities  of  various  types.  Most  of 
them  are  obvious  in  context,  and  there  is  no  particular  reason  to  go  through  them  now.  The  table  is  simply 
provided  as  an  easy-to-locate  reference  for  when  a  quick  review  of  notation  is  needed. 


Table  1  —  Notational  Conventions 


discrete-space  function 
particular  sample  of  above 
above  Fourier  transformed  on  n  index 
or  Fourier  transformed  on  n  and  t  both 

matrix 
three-vector 
two-vector 
three-vector  unit  vector 
real  or  real  vector 
integer  or  integer  two- vector 
constants  2.71828  . . .  and  \J  —  1 
column  vector  to  do  with  position 
row  vector  to  do  with  beamspace 


s°if) 

s°(/,f) 

B 

k,  x 
k,f,  n 

U,  V 

/,  w,x 

m,n,m 

e,  j 

t,  m,  n,  x 

k 


(the  function  as  a  whole;  the  “0”  could  be  any  label) 
(vector  n  could  be  any  allowed  index) 

(seldom  used) 

(uppercase  bold) 

(lowercase  bold) 

(lowercase  bold  italic) 

(pointy  hat  because  it  points) 

(not  from  center  of  alphabet) 

(from  center  of  alphabet) 

(the  usual  letters  but  in  a  roman  font) 

(especially  these  letters) 

(especially  these  letters) 


The  conventions  for  position-related  quantities  are  much  as  expected.  Variable  x  is  a  continuously  valued 
three-vector  position,  and  variable  n  is  a  discrete-space  sample  index  (see  next  section).  But  the  convention 
here  for  beamspace  quantities,  the  last  line  of  Table[l]  is  a  touch  tricky,  as  it  violates  the  Table [T] convention 
of  using  letters  from  the  center  of  the  alphabet  only  for  integer  valued  quantities.  In  particular,  this  report’s 
beamspace  is  closely  related  to  that  property  of  a  propagating  plane  wave  that  in  the  electromagnetics  lit¬ 
erature  is  typically  denoted  by  k  and  termed  its  wavenumber  vector,  wave  vector,  or  propagation  vector. 
For  this  variations  of  the  letter  k  and  nearby  letters  of  the  alphabet  are  used  for  most  beamspace  quantities. 
In  Section  3.3  below,  we  see  that  the  vector  k  itself  in  this  report  represents  the  same  information  as  that 
conventional  wavenumber  vector  but  in  a  slightly  different  way.  Furthermore,  our  use  of  the  discrete  Fourier 
transform  (DFT)  results  in  its  frequency  domain  being  a  subset  of  beamspace,  so  the  usual  DSP  convention 
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of  indexing  DFT  frequency  bins  with  a  vector  k  fits  this  report’s  notation  scheme  naturally  enough.  Finally, 
since  beamspace  is  actually  a  space  of  vector  spatial  frequencies,  our  spatial-frequency  vector  f  falls  into 
this  group  as  well. 


2.  LATTICE  BASICS 

2.1  The  Array  Plane:  Bases  and  their  Duals 

This  one  section  develops  the  only  linear  algebra  we  need  that  goes  beyond  ordinary  matrix  algebra. 
The  aim  is  to  go  slowly  and  carefully  here  but  to  depend  as  much  on  a  basic  visual  understanding  of  the 
geometry  of  the  plane  as  on  anything  very  formal. 

Certain  important  basics  about  the  geometry  of  planes  are  quite  familiar  and  easily  visualized.  Let’s 
think  about  the  array  plane  in  particular,  and  let’s  suppose  it  passes  through  the  origin.  For  now  we  need 
not  constrain  its  orientation.  Three  noncollinear  points  determine  a  plane,  so  take  the  origin  as  one  of  the 
three  points  and  choose  any  pair  of  nonzero,  noncollinear  three-vectors  b1  and  b2  in  the  array  plane  as  the 
other  two  points.  Figure  [T]  shows  an  example,  with  the  origin  in  the  center  and  with  the  half  arrowhead 
extending  clockwise  or  counterclockwise  from  a  vector’s  shaft  according  as  it  depicts  b1  or  b2,  conventions 
used  throughout  this  document.  Making  these  vectors  noncollinear  makes  them  linearly  independent,  which 
in  essence  means  that  any  point  p  in  the  array  plane  can  be  reached  by  summing  some  amount  a  of  vector 
b1  and  some  amount  (3  of  vector  b2  so  that  p  =  ab 1  +  /3b2.  Such  vectors  b1  and  b2  arc  basis  vectors  for  the 
array  plane,  and  the  pair  (b1,  b2)  of  two  vectors  is  a  basis  for  the  array  plane. 


-o.stF+o.gb2 


Fig.  1  —  Any  point  in  the  array  plane  can  be  written  as  a 
weighted  sum  of  basis  vectors  b1  and  b”,  shown  with  half 
arrowheads  extending  clockwise  and  counterclockwise  from  the 
shaft  respectively. 


The  whole  point  of  working  with  a  basis  in  this  setting  is  to  be  able  to  conveniently  specify  a  position 
like  p  in  the  plane  with  two  numbers  a  and  (3 ,  a  sort  of  coordinate  pair,  instead  of  having  to  specify  x,  y,  and 
z  components  of  p.  Also,  using  a  basis  means  we  need  not  be  concerned  with  which  combinations  of  x,  y, 
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and  2  lie  in  the  plane.  We  can  make  a  and  [3  anything  we  like,  and  p  is  in  the  array  plane  automatically.  We 
need  not  restrict  basis  vectors  b1  and  b2  to  have  unit  length  or  to  be  orthogonal  to  each  other.  While  such 


orthonormal  bases  are  quite  common  in  other  settings,  we  see  below  in  Section  2.2  that  here  we  have  good 
reasons  to  avoid  requiring  orthonormality.  It  is  enough  for  us  that  the  basis  vectors  are  in  the  array  plane 
and  are  linearly  independent,  i.e.,  not  collinear,  that  neither  is  just  a  scalar  multiple  of  the  other. 


Matrix  notation  is  very  helpful  in  working  with  basis  vectors.  When  the  vectors  are  all  column  vectors, 
the  sum  p  =  ab2+  /3b2  can  be  put  in  the  form 


column 
two-vector 
of  weights 


3x2 

basis  matrix 


When  the  vectors  are  instead  all  row  vectors, 


P 


row 

two-vector 
of  weights 


2x3 

basis  matrix 


Both  column  and  row  forms  are  used  heavily  below,  columns  for  position  vectors  and  rows  for  beamspace 
vectors.  Using  both  in  this  way  saves  us  many  transposes,  because  the  frequency-time  product  ft  that  occurs 
so  often  in  ordinary  Fourier  mathematics  has  its  counterpart  here  in  the  inner  product  of  a  beamspace  vector 
and  a  position  vector.  Here  product  (beamspace  vector)  (position  vector)  can  be  written  without  transposing 
either  quantity. 

Let  us  focus  first  on  column  vectors.  The  matrix-product  representation  of  p  of  ([[])  can  be  written 
p  =  Bin,  where  3x2  basis  matrix  B  has  basis  vectors  b1  and  b2  as  its  columns  and  where  real  column 
two-vector  w  contains  the  particular  weights  needed  to  produce  p.  But  given  B  and  p,  can  we  solve  for 
wl  We  certainly  cannot  just  write  w  =  B  1  p,  because  the  basis  matrix  B,  being  nonsquare,  cannot  be 
inverted!  But  in  fact  the  Gram  matrix  of  all  the  dot  products  of  the  basis  vectors,  given  by  B 1 B,  always  can 
be  inverted.  Before  we  try  to  use  this  fact,  let  us  see  why  it  is  reasonable  to  believe  it  is  true. 

Suppose  real  two-vector  u  is  given.  Must  there  exist  a  real  array-plane  three  vector  x  such  that  u  = 
Bt  x?  If  we  express  the  basis  vectors  as  lengths  times  unit  vectors,  as  b1  =  Ijb1]!  61  and  b2  =  ||b2||  b2,  the 
condition  u  =  B 1  x  becomes 


=  b1  •  x  or 

„  equivalently 

M2  =  bz  •  x  just 


7/1/||b1||  =  bl  -x 
M2 /  ||b2 1 1  =  b2  ■  x. 


The  quantities  on  the  far  right  are  just  the  components  of  x  in  the  b1  and  6 2  directions,  so  just  thinking 
visually,  the  answer  is  yes,  there  is  certainly  some  array-plane  x  that  extends  to  a  prescribed  degree  along 
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the  directions  of  the  basis  vectors,  and  it  is  in  fact  unique.  (The  linear  independence  of  the  basis  vectors 
is  critical!)  And  of  course  that  x,  simply  because  it  is  in  the  array  plane,  can  itself  be  written  as  x  =  Bin 
for  some  real  two-vector  w  of  weights  (even  though  we  don’t  quite  know  yet  how  to  actually  solve  for  w) 
so  that  u  =  B  1  x  =  B  1  But.  We  have  in  effect  argued  that  for  any  u  it’s  possible  to  invert  the  system 
u  =  B  1  Bm  of  lineal-  equations  to  find  unique  solution  w.  We  can  also  transpose  that  system  of  linear 
equations  and  write  it  as  m1  =  w  l  B  1  B.  Since  we  can  obtain  w  from  u  in  either  case,  we  can  effectively 
undo  a  multiplication  by  B 1 B  irrespective  of  whether  that  multiplication  was  done  on  the  left  or  right.  This 
is  the  very  meaning  of  the  statement  that  matrix  B 1 B  is  invertible. 

So  back  to  our  original  problem,  that  of  solving  p  =  Bin  for  in.  We  can  premultiply  everything 
by  B1  to  obtain  Brp  =  B  1  Bin,  and  now  we  can  use  the  inverse  that  we  know  must  exist  and  write 
(B 1  B)_1B 1  p  =  in.  Comparing  to  our  original  relationship  p  =  Bin  shows  that  (B1  B)_1B 1  inverts  left 
multiplication  by  B.  It  does  not,  however,  invert  right  multiplication  by  a  tall  B,  so  we  cannot  properly  term 
it  the  inverse  of  B.  A  different  terminology  is  called  for. 

The  Moore-Penrose  pseudoinverse  of  any  matrix  B,  which  is  calculated  in  Matlab  by  the  pinv  (  ) 
function,  is  defined  for  a  “tall”  matrix  B  like  ours  bf] 

B+A  (BtB)~1Bt.  (2) 

Now  if  we  had  three  basis  vectors  in  B,  the  largest  number  of  linearly  independent  vectors  one  can  ever 
have  in  threespace,  basis  matrix  B  would  be  square.  Definition  (|2])  would  simplify  to  yield  B+  =  B  ,  a 
true  inverse.  But  our  basis  matrix  B  has  only  two  basis-vector  columns,  so  that  simplification  is  not  possible 
here.  If  does  not  quite  always  invert  then,  what  exactly  does  it  do?  How  does  it  behave? 

Several  properties  of  the  pseudoinverse  are  important  to  us  here.  The  rows  of  BTon  the  right  in  definition 
([2])  are  just  b1T  and  b2T,  the  transposes  of  our  original  array-plane  basis  vectors.  If  we  for  the  moment 
write  the  inverse  of  the  Gram  matrix  as  D  =  (B  1  Bp 1  so  that  pseudoinverse  definition  (|2])  becomes  just 
B+  =  DB  r,  the  rows  of  B+  are  then  given  by 

(B+)r0wl  =  O1)1b1T  +  O1,2b2T, 

(B+)row  2  =  02, t  b1T  +  D2;2  b2T. 

The  four  elements  of  the  2x2  matrix  D  become  weights  used  to  combine  the  original  two  basis  vectors, 
now  oriented  as  rows,  to  make  two  new  vectors,  the  rows  of  B+.  A  simple  conclusion  follows. 

First  useful  property:  the  rows  of  B  +  are  each  in  the  array  plane. 

Next  let’s  ask,  are  these  rows  linearly  independent?  Vectors  are  linearly  independent  if  their  linear 
combinations  are  unique.  For  example,  since  the  basis-vector  columns  of  basis  matrix  B  are  linearly  inde¬ 
pendent,  whenever  both  p  =  Brr  and  p  =  Bn  are  true  it  must  be  that  u  =  v  as  well,  because  there  can 

'We  won’t  actually  need  the  pseudoinverse  of  any  wide  matrices,  but  if  we  did  we  would  use  the  easily  remembered  fact  that 
B+t  =  Bt+,  that  transposing  and  taking  the  pseudoinverse  commute  and  so  can  be  done  in  either  order.  For  a  wide  B  then,  we 
can  transpose  both  sides  of  this  to  obtain  B+  =  BT+T  and  then  apply  the  definition  of  pseudoinverse  to  tall  matrix  BT.  For  the 
record,  combining  the  steps  algebraically  yields  that  when  B  is  wide,  B+  =  BT(BBT)-1. 
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be  only  one  set  of  weights  that  can  be  used  to  combine  the  basis  vectors  to  reach  a  given  point.  Of  course, 
an  equivalent  way  to  say  this  is  that  if  B(rr  —  v)  =  0,  it  follows  that  u  —  v  =  0.  The  classic  mathematical 
definition  of  linear  independence  is  along  those  lines:  the  columns  of  B  are  linearly  independent  if  and  only 
if  Btu  =  0  implies  w  =  0. 

Applying  a  row  version  of  the  definition  of  linear  independence  to  the  pseudoinverse  then,  rows  of  B+ 
are  linearly  independent  if  and  only  if  /'B  =  0  implies  row  two-vector  f  =  0.  So  suppose  /B+  =  0.  By 
definition  Q  then,  (/(B1B)_1)Bi  =  0.  But  the  rows  ofB1  are  just  the  columns  of  B  and  so  are  linearly 
independent,  so  the  weight  vector  /(B 1 B)-1  =  0.  Right  multiply  by  Gram  matrix  B 1 B  and  we  have  that 
/  =  0,  revealing  another  fact. 

Second  useful  property:  the  rows  of  B+  are  linearly  independent. 

Another  useful  property  then  follows  from  the  first  two. 

Third  useful  property:  the  rows  of  B+ form  a  basis  for  the  array  plane. 

In  general,  the  rows  of  B+are  a  different  basis  than  the  one  we  began  with,  the  columns  of  B.  (Actually  they 
are  the  same  basis  vectors  if  but  only  if  b1  and  b1  are  orthonormal.)  The  rows  of  B  1  arc  just  as  legitimate  a 
basis  for  the  array  plane,  however.  A  simple  identity  is  the  key  to  the  usefulness  of  using  B+as  a  basis. 

Fourth  useful  property: 

B+B  =  I.  (3) 


Substitution  of  definition  ([2])  gives  the  proof  immediately:  B+B  =  (B 1 B)  1 B  rB  =  I. 

Our  basis-choice  strategy  will  be  to  use  B  and  B+ respectively  as  bases  for  array-plane  row  vectors  and 
column  vectors,  for  two  reasons.  First ,  identity  (|3])  then  always  makes  solving  for  the  weights  easy.  Multiply 
column  vector  p  =  Bn;  on  the  left  by  B+  and  use  the  identity  to  write  B+p  =  B  Bw  =  w.  Or  multiply 
row  vector  k  =  /B+  on  the  right  by  B  and  use  the  identity  to  write  kB  =  /B+B  =  /.  Second,  identity  ([3]) 
greatly  simplifies  products  of  the  general  form  (beamspace  vector)  (position  vector),  where  both  vectors  are 
in  the  array  plane.  Any  such  row  vector  (beamspace  vector)  can  be  written  as  / B+  for  some  weight  vector 
/,  and  any  such  column  vector  (position  vector)  can  be  written  as  Bn  for  some  weight  vector  n.  (Beginning 
in  Section [272] the  latter  weights  are  integers,  which  is  why  the  letter  n  seems  a  good  choice  here.)  Then  the 
product  simplifies  because  identity  ([3])  gives  (beamspace  vector)  (position  vector)  =  (/B+)(Bro)  =  fn. 
Our  two  bases  simply  drop  out,  an  enormous  convenience  mathematically.  In  honor  of  this  very  special 
property,  we  can  reasonably  call  the  columns  of  B  and  the  rows  of  B + dual  bases.  Each  basis  is  the  dual  of 
the  other.  (Some  authors  instead  term  them  reciprocal  bases.) 

Some  examples  of  bases  and  their  duals  are  shown  in  Fig.  [2]  along  with  several  properties  helpful  in 
sketching  dual  bases  without  actually  computing  a  pseudoinverse,  particularly  for  the  case  common  in  array 
work  in  which  basis  vectors  have  equal  lengths.  Proofs  of  those  properties  are  omitted  here  because  they  are 
tedious,  but  none  is  particularly  difficult.  Of  special  note  is  the  last  property  listed:  basis  vectors  of  equal 
lengths  result  in  equal-length  vectors  in  the  dual  basis  as  well,  and  the  stems  of  the  original  and  dual  “tees’,’ 
drawn  lightly  in  the  figure,  have  lengths  that  multiply  to  exactly  1/2.  Of  note  as  a  result: 
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Fig.  2  —  Properties  that  relate  array- 
plane  bases  (sketched  on  the  left)  and 
their  corresponding  dual  bases  (sketched 
on  the  right).  The  most  practically  useful 
of  these  properties  is  the  last: 

•  Always:  B+B  =  I.  Illustrated  in  the 
top  row. 

•  Always:  Spacing  basis  vectors  by  <j> 
spaces  dual-basis  vectors  by  comple¬ 
mentary  angle  180°—  <j>. 

•  Always:  Drawing  “tees”  with  tops 
connecting  vector  tips  and  with  stems 
from  the  origin  to  top  midpoints  yields 
a  basis  stem/top  orthogonal  to  the  dual¬ 
basis  top/stem.  Illustrated  in  the  sec¬ 
ond,  third,  and  fourth  rows. 

•  When  basis  vectors  are  orthogonal: 
lengths  are  reciprocal  to  the  lengths  of 
the  corresponding  dual-basis  vectors. 
Illustrated  in  the  second  row. 

•  When  basis  vectors  have  identical 
lengths:  dual-basis  vectors  also  have 
identical  (to  each  other)  lengths,  paral¬ 
lel  stems,  parallel  tops,  and  a  product 
of  stem  lengths  of  exactly  1/2.  Illus¬ 
trated  in  the  third  and  fourth  rows. 
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•  Once  orientations  of  equal-length  vectors  are  fixed,  vectors  in  B  and  B+have  lengths  inversely  pro¬ 
portional  to  each  other.  Halving  one  doubles  the  other,  and  so  forth. 

•  That  number  1/2  is  dimensionless,  so  the  units  of  the  vectors  in  B  and  B+must  be  inverses.  In  array 
work,  these  units  are  often  length  and  inverse  length  respectively,  representing  position  and — it  turns 
out — spatial  frequency,  or  those  units  can  be  normalized  out  so  that  all  vectors  arc  dimensionless.  (We 
follow  the  latter  approach  below.) 


Of  course,  identity  ([3])  gives  B+one  of  the  key  properties  of  an  inverse.  If  it  were  a  true  inverse,  however, 
we’d  also  have  BB+  =  I,  but  in  fact  this  cannot  be  true  when  B  is  nonsquare.  We  have  to  settle  for  a  weaker 
statement. 

Fifth  useful  property:  BB+  acts  like  an  identity  matrix  when  multiplying  array-plane  vectors. 

To  see  this,  suppose  column  three-vector  x  is  in  the  array  plane  and  write  it  in  terms  of  the  basis  and  a 
real  column  two- vector  w  of  weights  as  x  =  B  w .  We  know  this  can  be  done,  and  in  fact  we  know  now 
how  we  could  even  solve  for  the  weights,  though  we  need  not  do  so  here.  Instead,  just  substitute  into  the 
left  side  of  what  we  want  to  prove  and  watch  the  right  side  fall  out  through  the  use  of  identity  ([3])  above: 
(BB+ 

)X  “  (BB+)(Bu>)  =  B(B+B)m  =  B  w  =  x.  It  is  very  similar  if  we  start  with  a  row  three-vector  k 
expressed  in  the  dual  basis  using  a  two- vector  /  of  weights:  k  =  f  B+.  Now  k(BB+)  =  (/B+)(BB+)  = 
/(B+B)B+  =  /B+  =  k.  It  would  be  very  reasonable  to  view  B+Bx  =  x  and  kB+B  =  k  as  dual 
properties. 

2.2  Array-Plane  Lattices  and  their  Duals 

In  the  discussion  above  an  arbitrary  array- plane  column- vector  position  p  often  takes  form  B  w  using 
a  real  two-vector  w  of  weights  and  a  3  x  2  basis  matrix  B  that  has  basis  vectors  as  columns.  Now  let 
us  consider  a  very  special  case,  that  of  integer  weights.  The  set  of  all  possible  integer-weighted  linear 
combinations  of  a  set  of  basis  vectors  is  a  point  lattice  or  just  a  lattice.  We  write  BZ2  for  the  lattice 
comprising  column  vectors  of  the  form  Bn  with  n  an  arbitrary  column  two-vector  of  integers — recall  that 
we  write  just  n  G  to  specify  such  n  concisely.  We  write  Z2B  '  for  the  lattice  comprising  row  vectors 
of  the  form  fcB  1  where  k  e  Z2,  that  is,  where  k  is  an  arbitrary  row  two-vector  of  integers.  Let  us  term 
integer  column  two-vector  n  and  integer  row  two- vector  k  the  index  vectors  that  specify  lattice  points  B  n 
and  A:B  2  The  lattices  BTrand  Z2B+ constructed  from  dual  bases  are  dual  lattices. 

The  left  and  right  columns  of  Fig.  [3]  show  basis  vectors  and  points  of  lattice  BZ2  and  its  dual  Z2B+ 
respectively.  Matrix  B  is  the  same  for  both  top  plots,  so  the  bases  shown  in  those  plots  are  duals.  A  different 
matrix  B  was  chosen  for  the  bottom  plots,  but  it  is  the  same  for  both  of  them,  making  the  bases  shown 
in  the  bottom  row  duals  as  well.  The  lattices  in  the  top  and  bottom  rows  are  the  same,  even  though  the 
index  vectors  for  corresponding  points  are  not,  so  it  is  clear  that  a  lattice  basis  is  not  unique.  In  fact,  we 
could  correctly  claim  that  the  lattices  in  the  upper-left  plot  and  the  lower-right  plot  are  duals  even  though 
the  bases  shown  and  used  to  construct  them  are  not.  Dual  bases  generate  dual  lattices,  but  dual  lattices  can 
be  constructed  from  bases  that  are  not  dual. 

The  dual  lattices  shown  are  hexagonal  lattices,  lattices  on  equilateral-triangular  grids.  Hexagonal  lattices 
with  the  nearest-neighbor  distances  shown,  1  j\f?>  on  the  left  and  2  on  the  right,  arc  classic  in  array  work. 
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Fig.  3  —  Left  column:  a  hexagonal  lattice  B$2  with  l/y/3  nearest-neighbor  spacing.  Right  column: 
the  associated  dual  lattice  Z12B+ is  also  hexagonal  but  with  a  nearest-neighbor  spacing  of  2.  Different 
basis  matrices  B  are  used  in  the  top  and  bottom  rows,  illustrating  that  bases  are  not  unique  and  that 
dual  lattices  can  be  generated  from  nondual  bases  as  well  as  from  dual  bases. 
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A  basis  of  equal-length  vectors  separated  by  60°  has  as  its  dual  a  basis  of  equal-length  vectors  separated  by 
120°,  but  these  basis-vector  lengths  produce  “tee”  stem  lengths  very  different  in  the  top  and  bottom  rows. 
The  top-row  lengths  of  1/2  and  1  are  easier  to  remember,  so  in  this  document  that  basis  geometry  is  favored. 

Given  a  basis  matrix  B,  the  easiest  way  to  check  whether  it  generates  one  of  the  hexagonal  lattices 
shown  is  to  compute  Gram  matrix  B1  B.  Its  diagonal  elements  are  the  squared  lengths  of  the  basis  vectors, 
and  given  that  those  are  the  same  for  the  two  vectors,  the  off-diagonal  elements  are  that  squared  length  times 
the  cosine  of  the  angle  between  the  vectors.  This  works  out  nicely  here  because  the  cosines  of  60°  and  120° 


are  1/2  and  —1/2  respectively.  The  four  bases  in  Fig.0then  should  have  these  four 

Gram  matrices. 

Fig-a 
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Rectangular  lattices  with  nearest-neighbor  distances  of  1/2  and  2  are  also  common  in  array  work  and 
have  bases  and  dual  bases  with  gram  matrices 


BtB  =  - 


1  O' 

0  lJ  ’ 


B+B+t  =  4 


1 

.0 


O' 

1. 


though  rectangular  lattices  are  not  be  considered  further  in  this  document. 

One  simple  consistency  check  on  the  above  matrices  is  provided  by  the  fact  that  the  Gram  matrix 
B  1  B  of  a  basis  matrix  B  and  the  Gram  matrix  B+B+T  of  dual  basis  matrix  B+  are  always  inverses 
because,  substituting  pseudoinverse  definition  ([2])  and  canceling  adjacent  inverse  matrices,  B+B+T  = 
((BiB)_1Bt)  (B(B1B)_1)  =  (BtB)-1.  For  simplicity  then,  it  is  natural  to  write  the  Gram  matrix 
of  dual  basis  matrix  B+as  just  (BTB)-1. 

2.2.1  Constructing  a  Basis  for  Array  Work 

Let  us  construct  the  matrix  B  with  column  basis  vectors  as  shown  in  the  upper-left  plot  of  Fig.  [3] 
We  could  pick  any  orientation  for  the  array  plane,  but  let  us  pick  one  that  makes  the  construction  simple. 
Visualize  three-space  unit  vectors  in  the  usual  x,  y,  and  z  directions  and  notice  that  if  we  connect  their  tips 
and  view  it  from  a  point  far  out  in  the  [1  1  1]  direction,  we  see  an  equilateral  triangle.  We  can  therefore 
obtain  a  pair  of  vectors  at  the  desired  60° angle  by  subtracting  unit  vectors,  following  which  we  can  normalize 
to  get  the  desired  lengths.  These  steps  yield 

make  make 

desired  unit 


(4) 
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We  can  immediately  assemble  basis  vectors  b1  and  b2  into  basis  matrix  B  and  compute  dual-basis  matrix 
B+  from  it: 


B  = 


1 


'-1 

0 

o' 

-1 

,  B+=  v^73 

■-2  i  r 
.1-2  1. 

1 

1 

(5) 


It  is  easy  to  check  that  B+B  =  I,  that  basis  matrix  B  has  Gram  matrix  B 1 B  =  |  [ 2  2  ] ,  and  that  dual-basis 
matrix  B+has  Gram  matrix  B+B+T  =  2  [  2  ~^]  =  (B  1  B)~  1  as  required. 


2.2.2  Plotting  Axes:  Azimuth  and  Elevation 


The  boresight  vector  is  normal  to  the  array  plane  and  generally  presumed  in  array  discussions  to  be 
towards  the  horizon  or  be  at  some  modest  array  “tilt  back”  angle  above  the  horizon.  Most  often  we  view  the 
array  plane  in  diagrams  from  behind,  looking  in  the  boresight  direction. 


In  construction  (J4])  the  geometry  is  viewed  from  far  off  in  the  [1  1  1]  direction,  so  we  are  actually 
looking  in  the  [—1  —1  —  l]  direction.  We  can  see  in  Q  that  all  our  basis  vectors,  the  rows  of  B  and  the 
columns  of  B+,  are  orthogonal  to — have  zero  dot  products  with — this  vector,  so  here  let  us  take  the  direction 
of  [— 1  —1  — l]  as  the  boresight  direction. 


The  rightward  and  upward  directions  on  the  page  in  construction  Q  are  then  named  the  azimuth  and 
elevation  directions  simply  because  if  we  start  with  a  boresight  vector  extending  towards  the  horizon  and 
rotate  it  towards  vectors  in  those  azimuth  or  elevation  directions,  we  respectively  increase  our  boresight 
vector’s  compass  direction  or  angle  above  the  horizon.  We  can  formalize  this  by  constructing  unit  vectors 
in  the  three  orthogonal  directions,  azimuth,  elevation,  and  boresight.  In  Q  the  azimuth  direction  appears  to 
be  the  z  direction,  but  if  we  start  with  an  arbitrary  vector  in  that  direction,  for  example  [0  0  3] ,  it  is  not 
orthogonal  to  the  boresight  direction  as  required.  We  can  correct  this  by  adding  a  boresight  component,  in 
this  case  just  [—1  —1  — l] ,  so  that  the  total,  [—1  —1  2] ,  is  indeed  orthogonal  to  boresight.  Similarly, 

the  elevation  direction  appears  in  ([4])  to  be  in  the  direction  of  [l  —  1  0],  with  equal  components  in  the 

x  and  —y  directions,  and  we  find  this  vector  is  orthogonal  to  boresight  without  adjustment.  Let  us  now 
construct  a  matrix — call  it  P  for  plotting — with  our  three  orthogonal  vectors  as  columns  and  with  adjusting 
scale  factors  included  as  a  diagonal  post-multiplying  matrix  to  give  each  column  unit  length. 


in  azimuth  direction 


I  I 


in  elevation  direction 


r-i  i  -ii 

1 - 

o 

o 

~f> 

1 _ 

P  = 

-l  -i  -i 

o 

H“> 

o 

2  0  -lj 

s 

L  0  0 

normalization  to  unit  length 


It  is  easy  to  verify  that  PTP  =  PPT  =  I,  making  it  an  orthogonal  matrix  as  intended.  Now  array- 
plane  points  like  column  vector  x  =  Bn  or  row  vector  k=/B+can  be  projected  onto  the  columns  of  P  to 
obtain  their  az-el-boresight  coordinates  as  P  1  x  =  P  1  Bn  or  kP  =  f  B+P  respectively.  If  we  have  done 
everything  correctly,  the  last  component,  boresight,  should  be  zero.  The  first  and  second  coordinates,  az  and 
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el,  are  plotted  on  horizontal  and  vertical  axes  to  make  array-plane  plots.  If  we  wish,  we  can  precompute  the 
matrix  products 
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Of  course  we  can  omit  the  boresight  column  from  P  if  we  wish  and  make  P  a  3  x  2  az-el  plotting  matrix  only, 
making  both  P 1 B  and  B+P  into  2x2  matrices.  When  working  with  the  hexagonal  lattice,  many  authors 
avoid  three-dimensional  representations  and  work  with  two  array-plane  coordinates  only,  effectively  using 
these  2x2  versions  of  P1  B  and  B  P,  which  are  (easily  verified  to  be)  inverses,  as  their  basis  and  dual-basis 
matrices.  At  first  glance  it  seems  attractive  to  work  in  threespace  for  the  naturalness  with  which  the  60°  and 
120°  angles  can  be  obtained  so  that  radicals  are  confined  to  scalar  factors  out  front,  rather  than  have  radicals 
scattered  across  various  elements  of  the  2x2  versions.  But  it  actually  makes  no  difference,  because  one 
plots  az  and  el  and  so  must  compute  azimuth  and  elevation  for  plotting  using  P  1  B  and  B+P  anyway,  and 
in  other  contexts  the  B  and  B+ matrices  most  often  simply  cancel. 


2.3  Sublattices  and  their  Duals 


What  happens  when  we  use  points  of  a  lattice  as  basis  vectors  to  generate  a  new  lattice? 

Suppose  we  begin  with  the  bases  B  and  B+  and  the  lattices  B^2  and  2rB  1  they  generate,  exactly  as 
pictured  in  the  top  row  of  Fig.  pi  and  choose  two  points  of  the  lattice  B^2  on  the  left  to  use  as 

.  The  upper  half  of  Fig^h  shows  one  such  choice,  the  specific  points  B  [  [  ]  and  B  [  .  We  can 

make  these  the  columns  of 


B 

1 

,  B 

'-r 

=  B 

1  -1 

=  BR  with  R  = 

'1  -1' 

.1. 

.  2. 

.1  2. 

.1  2. 

Choosing  points  of  the  original  lattice  as  gives  us  a  BR. 

The  lower-left  plot  in  Fig.  [d]  shows  the  BR|2  generated  by  this  BR.  A 

point  in  that  is  of  the  form  BRn  with  n  which  can  be  looked  at  in  two  ways.  It  is  of 

course  basis  matrix  BR  times  n  e  which  is  what  makes  it  a  point  of  BR|2,  but  it  is  also 
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Fig.  4  —  Top  left:  make  R  a  nonsingular  integer  resampling  matrix  to  make  ,  the 

columns  of  BR,  be  points  of  an  existing  lattice  Top  right:  the  original  dual  basis  vectors 

are  integer  combinations  of  the  ,  the  rows  of  R_1B+.  Bottom  left:  the 

generate  a  BRf  of  original  lattice  B^2.  Bottom  right:  the  dual  of  the 

is  a  S2R~1B+  of  the  original  dual  lattice  Z2B+.  Summary  of  the 

relationships: 

original  lattice  Bf?  <^Ual  Z2B+ 

U  density  ratio  |  R|  (T 

sublattice  of  original  lattice  BR|2  < - >  Z2B,~1B+ 
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basis  matrix  B  times  Rn  G  which  makes  it  simultaneously  point  of  the  original  lattice  B.?p2  In  this  way, 
every  point  in  BR|f2  is  also  in  B^2,  so  we  write  BR^2  C  B^2and  say  the  BR,|f2  is  a  sublattice 

of  the  original  lattice  B|f2.  This  only  works  because  R  is  a  resampling  matrix,  a  nonsingular  square  integer 
matrix.  (“Nonsingular”  here  corresponds  to  the  being  linearly  independent  Of  course 

in  Fig.  [4]  we  could  have  chosen  completely  different  points  for  our  in  a  way  that  yielded 

the  same  ,  so  resampling  matrices  are  not  unique.  Here  R  =  [ J  as  pictured,  but  BR  with 

R  =  [  “  J  2  ] ,  for  example,  would  have  generated  the  same  sublattice. 

What  does  the  dual  of  a  sublattice  look  like? 

The  (row)  basis  B+  and  the  lattice  Z2B  :  of  row  vectors  it  generates  are  duplicated  from  the  upper 
right  in  Fig.  [4]  to  the  upper  right  in  Fig.  [3j  where  the  dual  of  our  is  then  added.  That  dual 

is  easily  computed  using  a  simple  algebraic  property  of  the  pseudoinverse.  Here  the  derivation  of  that 
property  applies  pseudoinverse  definition  (|2])  to  BR  in  the  first  line,  distributes  transposes  and  inverses 
across  products  in  the  second  and  third  lines,  cancels  inverses  in  the  fourth  line,  and  recognizes  the  right 
side  of  definition  Q  at  the  end: 


(BR)+=  ((BR)tBR)_1(BR)t 
=  (RtBtBR)"1RtBt 
=  R“1(BtB)1R"tRtBt 
=  R“1(BtB)“1Bt 

=  R-1B+.  (6) 

The  similarity  in  form  to  the  matrix-algebra  identity  (UV)'1  =  V_1U_1  makes  this  easy  to  remember. 
The  basis  to  our  ,  the  columns  of  BR,  comprises  just  the  rows  of  R_1B+  or,  using  the 

Fig. [^labels,  [l  OjR^B+and  [0  ljR-^t 

Hidden  within  identity  ([6])  is  an  important  relationship.  Look  at  the  basis  vectors  of  the  original  dual 
lattice.  In  Fig.  [3J  the  upper  and  lower  of  these  basis  vectors  were  labeled  [0  l]B+and  [l  0]B+  because 
these  vectors  are  just  the  two  rows  of  B+.  But  since  RR_1  =  I,  we  are  free  to  insert  RR  1  into  these 
expressions  anywhere  that  matrix-vector  sizes  are  compatible.  We  can  say,  for  example,  that  the  two  basis 
vectors  are  [0  l]RR^1B+and  [l  0]  RR_1B+,  revealing  them  to  be  [0  l]R  =  [l  — l]  and  [l  0]R  = 
[l  2]  times  R^1B+,  which  is  how  they  are  labeled  in  Fig.  |d]  So  in  addition  to  being  basis  vectors  of  our 
original  dual  lattice,  they  are  also  points  in  some  Z2R_1B+ generated  by  R-1B+, 

the  pseudoinverse  and  therefore  dual  of  BR.  We  have  discovered  that  Z2B+  C 

Z2R"1B+. 

This  Z2R_1B+  is  shown  on  the  lower  right  in  Fig.  [4]  On  the  left,  before  we  took  duals,  the 

was  a  sublattice  of  the  old,  but  on  the  right,  where  we  are  looking  at  their  duals,  this  is  reversed. 
The  old  lattice  is  a  sublattice  of  the  ,  or  the  is  a  superlattice  of  the  old.  The  dual  and  sublattice 

relationships  are  summarized  at  the  end  of  the  Fig.[4]caption. 

2We  know  that  BRw  =  0  if  and  only  if  Rw  =  0,  because  the  columns  of  B  are  linearly  independent,  and  if  we  choose  linearly 
independent  ,  we  also  know  that  BRw  =  0  if  and  only  if  w  =  0.  Combine  the  two  statements  and  we  have  that 

Rw  =  0  if  and  only  if  w  =  0,  which  for  a  square  R  is  just  the  statement  that  R  is  nonsingular. 
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Figure  [4]  also  illustrates  a  convention  used  throughout  this  document.  On  the  left  a  lattice  is  always 
plotted  before  its  sublattice.  As  each  is  plotted,  its  dual  is  plotted  on  the  right,  where  a  lattice  is  therefore 
plotted  before  its  superlattice.  Larger  dots  are  always  used  for  less-dense  lattices,  so  on  the  left  larger 
sublattice  dots  obscure  some  points  of  the  original  lattice,  forcing  us  to  remember  that  the  obscured  points 
are  there.  There  is  no  such  difficulty  on  the  right. 

2.3. 1  The  Density  of  a  Lattice  in  the  Plane 

How  do  we  determine  the  density  of  a  lattice,  i.e.,  the  number  of  its  points  per  unit  area  in  the  plane? 

Consider  the  example  lattice  on  the  left  in  Fig.  [5]  where  lattice  points  are  shown  as  enormous  dots  to 
simplify  this  discussion  and  where  crosshatching  through  the  lattice  points  has  tiled  the  plane  with  paral¬ 
lelograms.  The  area  of  one  of  these  parallelogram  tiles  is  the  fundamental  volume  of  the  lattice  (where  the 
term  “volume”  refers  to  the  possible  generalization  to  more  than  two  dimensions).  The  four  pieces  of  a  dot 
in  each  parallelogram  total  to  exactly  one  dot,  so  we  can  regard  the  tiles  as  unit  cells  and  take  the  density  of 
the  lattice  as  the  reciprocal  of  its  fundamental  volume. 


Fig.  5  —  Left:  division  of  the  plane  into  unit  cells  of  the  lattice. 

Right:  The  derivation  of  the  fundamental  volume  of  lattice  B|Fas  jECBj1/2. 


The  basis  vectors  are  sketched  Fig.  [5]  with  less  than  90°  between  them,  a  property  we  depend  on  here  in 
deriving  the  fundamental  volume.  If  the  actual  angle  is  more  than  90°,  just  replace  one  of  the  basis  vectors 
with  its  negative  to  make  the  angle  between  them  less  than  90°.  It  becomes  clear  below  that  we  can  actually 
just  forget  about  having  done  this  step,  because  the  final  algebraic  result  is  not  affected  by  negating  either 
basis  vector. 

To  obtain  the  fundamental  volume,  consider  the  enlarged  unit  cell  on  the  right.  It  is  a  parallelogram 
comprising  a  square  and  two  triangles.  If  the  unshaded  triangle  is  clipped  off,  slid  over,  and  joined  to  the 
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shaded  triangle  hypotenuse-to-hypotenuse,  the  original  square  and  two  triangles  become  one  large  rectangle. 
In  the  sketch  the  long  side  of  the  rectangle  has  length  1 1  b 1 1 1 .  The  length  of  the  short  side  of  the  rectangle  is 
the  same  as  that  of  the  unlabeled  side  of  the  shaded  right  triangle  shown.  The  hypotenuse  of  that  triangle 
has  length  |  j  b  2 1 1 ,  and  the  short  side  has  length  that  is  j  ust  the  dot  product  of  b 2  with  a  unit  vector  b 1  in 

the  direction  of  b1.  The  Pythagorean  theorem  gives  the  missing  triangle  side  length  as 


2 


so  the  large  rectangle  area,  the  area  of  the  unit  cell  or  fundamental  volume,  is 


l|b1||y//||b2||2  -  (^p^)2=  ((||b1||2||b2||2  —  (b1-  b2)2))  1/2  =  |det(BTB)|1/2. 

For  brevity  we  ordinarily  write  this  fundamental  volume  as  B'1  B  | 1  /  j  and  the  corresponding  lattice  density 
as  IB1  B|~1'/2.  The  determinant  is  not  affected  by  changing  the  sign  of  a  basis  vector,  so  this  is  a  universal 
result  when  the  basis  vectors  arc  columns  of  a  matrix  B. 


When  the  basis  vectors  arc  rows  of  a  matrix,  the  transpose  must  go  on  the  other  factor  to  create  the 
desired  dot  products  and  squared  lengths,  so  the  fundamental  volume  of  dual  lattice  Z2B+  is  jB+B+T|1/2. 
The  determinant  is  applied  to  the  Gram  matrix  of  dual  basis  B+,  and  of  course  that  Gram  matrix  can  be 
written  (B  1  B)  \  The  fundamental  volume  of  dual  lattice  Z2B  1  therefore  is  B  1  B  l,/2  and  this  dual’s 
density  is  |B  1  B i 12  The  lattice  and  its  dual  have  reciprocal  fundamental  volumes  and  reciprocal  densities, 
and  the  fundamental  volume  of  one  is  the  density  of  the  other. 

A  sublattice  BR^2  has  fundamental  volume 

|(BR)t(BR)|1/2  =  |RtBtBR|1/2  =  (|R|  |BtB|  |R|)1/2  =  |BtB|1/2|R|. 

This  is  an  integer  factor  |R|  larger  than  the  fundamental  volume  of  lattice  BZ2  so  sublattice  BR^2  has  a 
density  lower  than  that  of  lattice  B|Z  by  a  factor  of  |R|. 

Dual  lattices  have  reciprocal  fundamental  volumes,  so  lattice  Z2R_1B+,  a  superlattice  of  lattice  Z2B  * 
and  the  dual  of  lattice  BR|2,  has  fundamental  volume  B'B  i,/2/[R|,  a  factor  of  |R|  less  than  that  of 
lattice  Z2B+. 


3,  A  SIGNAL-PROCESSING  ENGINEER’S  VIEW  OF  THE  ELEMENT  SYSTEM 

Let’s  establish  the  array  plane  by  supposing  we  have  a  normalized,  dimensionless  3x2  basis  matrix 
B  such  that  any  physical  position  in  the  array  plane  can  be  represented  as  a  linear  combination  of  the 
basis  vector  columns  of  matrix  AB,  where  A  is  some  length-dimensioned  scaling  constant  (not  necessarily 
wavelength,  as  we  shall  see).  We  eventually  discover  rationales  below  for  choosing  both  B  and  A. 

We  can  now  develop  a  convenient  characterization  of  the  element  outputs. 
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3.1  A  Simple  and  Effectively  Periodic  Element  System 

As  a  thought  experiment,  suppose  there  is  nothing  in  the  universe  besides  incident  electromagnetic  fields 
(and  something  far,  far  away  that  creates  them)  and  a  receive  antenna  structure  driving  a  receiver  front  end 
that  in  turn  yields  a  single  output.  For  simplicity  suppose  the  structure  uses  only  ordinary  linear  materials 
and  that  the  electronics  arc  LTI,  i.e.,  just  linear  amplification  and  filtering.  Now  physically  translate — slide 
without  changing  orientation — this  entire  system  through  space  by  some  vector  x  and  write  the  output  signal 
as  a  function  of  that  translation  vector  as  y( t,  x).  Given  two  arbitrary  field  configurations  such  that 

(incident  field  configuration  1)  (t,  x)  yields  y1( t,  x), 

(incident  field  configuration  2)  (t,  x)  yields  y2( t,  x) , 

then  for  arbitrary  constants  a  and  6, 


a(field  configuration  1) (t,  x)  +  /?(field  configuration  2)(t,  x)  yields  ay1( t,  x)  +  (3y2(t,  x) 

just  by  the  lineality  of  Maxwell’s  equations  and  the  lineality  of  the  front  end,  so  the  system  is  linear.  In 
addition. 


(field  configuration  l)(t  —  t(  x)  yields  yl(t  —  t(  x). 

for  any  time  delay  t(  by  the  time  invariance  of  Maxwell’s  equations  and  the  front  end.  So  the  system  is  time 
invariant.  But  does 


(field  configuration  1 )  (t,  x  —  x')  yield  y1( t,  x  —  x7), 

for  arbitrary  vector  x'?  Is  the  system  space-invariant?  Indeed  it  is,  because  we  defined  the  spatial  argument 
of  y( t,  x)  as  the  vector  specifying  the  amount  by  which  we  slide  the  entire  structure  through  space.  If  we 
move  the  fields  east  by  one  meter  (by  moving  the  far-away  sources  east  one  meter)  and  move  the  antenna 
structure  east  by  one  meter  also,  the  output  function  of  time  is  unchanged. 

If  we  want  to  use  a  linear,  time-invariant,  space-invariant  (LTSI)  antenna  and  front-end  system  in  a  way 
that  takes  advantage  of  space  invariance,  what  are  we  to  do?  We  can’t  very  well  be  constantly  moving  this 
antenna  system  all  about,  can  we? 

We  almost  can,  in  fact.  What  we  can  do  is  create  what  we  pretend  is  a  two-stage  system.  The  first  stage 
is  lineai;  time-invariant,  and  space-invariant,  and  the  second  stage  comprises  spatial  sampling  at  the  points 
of  lattice  ABE2  Now  we  only  have  to  be  able  to  translate  the  first-stage  structure  by  an  arbitrary  x  €  ABE2, 
because  the  second-stage  sampler  never  asks  to  see  the  first-stage  output  for  other  translations.  The  one  way 
to  make  this  easy  to  do  is  to  make  the  structure  periodic.  We  make  it  invariant  to  translation  by  an  arbitrary 
x  €  ABE2  so  that  moving  it  changes  nothing.  We  are  free  to  pretend  we  moved  it,  as  no  one  can  ever  know 
we  didn’t.  It  is  like  the  moment  in  the  famous  mirror  scene  of  the  1933  Marx  Brothers'  film  Duck  Soup  in 
which  Haipo,  pretending  to  be  Groucho's  reflection  in  a  nonexistent  mirror,  is  caught  off  guard  by  Groucho 
suddenly  spinning  around  and  just  pretends  he  spun  as  well,  with  Groucho  none  the  wiser! 
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To  create  a  periodic  structure,  in  principle  we  just  begin  with  some  physical  antenna-element  structure 
that  tits  in  a  unit  cell  of  lattice  AB|i2.  Let’s  suppose  this  structural  cell  has  a  single  feedpoint  at  nominal 
location  x  =  0,  and  let’s  suppose  this  feedpoint  drives  LTI  front-end  amplification  and  filtering  that,  in 
principle  anyway,  we  can  suppose  is  also  contained  in  the  unit  cell.  For  convenience  let’s  assume  that  this 
front-end  processing  includes  enough  buffering  that  its  output  can  be  observed,  measured,  passed  on  to  other 
processing,  and  so  forth,  without  affecting  the  fields. 

Copying  this  unit-cell  structure  an  infinite  number  of  times,  offsetting  copies'  physical  positions  from 
the  original  by  amounts  A  Bn  with  n  £  rf?,  creates  a  periodic  structure  as  desired.  We  can  take  the  original, 
uncopied  front-end  output,  the  one  nominally  at  feedpoint  location  x  =  0,  as  the  output  y(t,  0)  of  the 
unmoved  first-stage  system.  Then  we  can  write  the  second-stage  output  as  (t)  4  y(  t,  A  Bn,)  for  sample 
index  n  For  simplicity  we  can  just  refer  to  this  as  the  output  of  the  element  “at”  physical  position  ABn 
or  the  element  indexed  by  vector  n  or  just  element  n.  Flere  the  superscript  zero  is  just  an  identifying  flag 
referring  to  these  processed  element  outputs.  Going  forward  ABn  and  XBf?  arc  respectively  termed  the 
physical  element  position  and  physical  element-position  lattice,  and  A-normalized  versions  Bn  and 
are  just  an  element  position  and  the  element-position  lattice  respectively. 

The  only  paid  of  our  thought  experiment  that  is  apparently  impractical  to  approximate  in  practice  is  its 
infinite  spatial  extent.  But  if  we  design  our  back-end  processing  to  use  only  a  finite  number  of  the  element 
outputs,  the  remaining  purpose  of  those  unit  cells  with  unused  outputs  is  to  provide  a  periodic  electro¬ 
magnetic  environment  for  the  unit  cells  with  outputs  we  actually  do  use.  This  is  important  enough  that  in 
general  we  should  not  simply  discard — fail  to  physically  construct — unit  cells  with  unused  outputs.  But  we 
can  provide  the  desired  electromagnetic  environment  for  the  outputs  used  by  keeping  nearby  unused  cells, 
the  antenna  portions  anyway,  and  replacing  the  missing  electronics  with  equivalent  termination  impedances. 
These  guard  elements  or  dummy  elements  make  our  array  seem  infinite,  at  least  as  far  as  we  can  tell  by 
inspecting  the  finite  number  of  outputs  actually  available  to  us. 

Finally,  let’s  use  the  term  elements  from  here  forward  for  the  combination  of  the  physical  antenna  ele¬ 
ments  and  any  LTI  processing  that  is  either  actually  performed  identically  on  all  physical-element  outputs  or 
that  can  be  referred  to  the  feedpoints  and  incorporated  into  our  notion  of  the  elements  for  analysis  purposes. 
This  would  naturally  include  any  frequency  responses  associated  with  per-element  front-end  processing,  but 
it  can  also  include  frequency  responses  associated  with  frequency  conversion  and/or  DSP  that  is  actually 
performed  downstream,  after  front-end  processing,  as  long  as  there  is  an  equivalent  front-end  operation  that 
can  be  included  in  the  element  model.  This  referral  is  more  than  a  minor  convenience,  as  the  simplified 
approach  to  analysis  taken  below  requires  that  certain  frequency  responses  be  included  in  the  elements  for 
proper  array-system  functioning  to  be  easily  seen. 

3.2  Fourier  Representation  of  the  Element  Outputs 

3. 2. 1  2D  Discrete-Space  Fourier  Transforms 

Given  some  function  s(ni,ri2)  of  integer  indices  n\  and  712,  one  certainly  might  treat  ri2  as  a  fixed 
parameter  and  take  a  discrete-time — here  it  actually  represents  space — Fourier  transform  on  n\  to  obtain,  if 
we  use  f\  as  the  frequency  variable, 


OO 

S(/i,n2)=  ^2  s(m,n2)e~j27rfini  . 

m=— 00 
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There  is  nothing  to  prevent  us  from  then  proceeding  to  take  f\  as  a  fixed  parameter  and  further  transforming 
s(/i,  712)  on  ri2  using  a  new  frequency  variable  /2,  resulting  in,  if  we  leave  the  infinite  limits  as  understood, 

S(h,  f2)  =  E  s(/i>  «2)  E  ®("1.  "2)  e-j27r(/mi+/2n2)  • 

712  n  2  ni 


Except  for  the  order  of  summation,  which  is  irrelevant  given  practical  restrictions  on  the  function  s(m,  712), 
the  same  result  exactly  would  be  obtained  if  we  did  the  transforms  in  the  other  order.  So  we  may  as  well 
define  vectors  n  =  [^  ]  and  /  =  [/j  f2 ]  and  write  the  whole  thing  using  a  vector  notation  as 

S(f)  =  Y,sne~i2nfn-  (7) 


In  the  interests  of  uncluttered  brevity,  in  this  document  sums  over  integer  column  vectors  or  integer  row 
vectors  arc  by  default  over  all  elements  of or  1?  respectively,  so  if  n  is  used  in  the  summand  as  a  column 
two-vector  of  integers, 

X]  means  ^  ,  which  is  just  ^  ^  °r  ^  ^  with  n  \n2  ] 

n  ne|2  m  n2  n2  m 


and  if,  as  in  later  sections  below,  k  is  used  in  the  summand  as  a  row  two-vector  of  integers, 


means  ,  which  is  just  or 

k  fees2  fci  k2 


EE 

k2  fci 


with  k  =  [fci  ,k2\. 


These  sums  arc  formally  infinite  but  in  reality  finite,  because  only  finitely  many  of  their  terms  arc  ever 
nonzero.  There  arc  no  convergence  questions  with  finite  sums,  so  their  terms  can  be  added  in  any  order,  and, 
consequently,  their  component  sums  on  individual  vector  components  can  be  ordered  in  any  way  desired. 

The  inverse  transforms  arc  handled  similarly,  except  that  integrals  replace  sums.  Inverse  transforms  of 
S(fi,  /%)  on  frequency  variables  f\  and  fi  can  be  done  in  either  order.  Inverse  transforming  on  /2  yields 


s(/i,n2)=  [  S(hJ2)e^f™dh 

J  period2 

where  region  of  integration  “period2”  is  any  one  period  of  the  integrand,  for  example  a  unit  interval  like 
—  1/2  to  1/2  or  0  to  1.  Then  inverse  transforming  on  f\in  an  exactly  similar  way  yields 

s(rai,n  2)=  [  S  {h,n2)^^dh  =  f  [  S(hJ2)e^^+^  df2dff. 

J  period  J  period  1  J  period2 

Again  we  may  as  well  use  a  vector  notation  and  write 

sn=  [  S(f)e^ndf, 

J  □ 


(8) 
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where  we  adopt  the  convention  that  df  means  differential  area  d/2  dfi  and  take  the  integral  to  be  a  double 
integral  because  two  variables  are  actually  integrated.  Writing  □  for  the  region  of  integration  here  is  intended 
as  shorthand  for  any  period  in  /.  This  could  be  (but  doesn’t  have  to  be)  just  a  rectangle  extending  over  a 
period  in  f\  and  a  period  in  /2.  Certainly  the  simplest  choice  is  a  unit  interval  in  each  variable,  which  of 
course  makes  a  square  in  the  (/j,  /2)  or  /  plane. 

Together  ([7])  and  ([8])  define  2D  discrete-space  Fourier  transform  pair  sn  ~  S(f). 

3.2.2  The  Element  Outputs  in  2D  and  3D  Frequency  Domains 

Inverse  discrete-space  Fourier  transform 

4(t)  =  [  S°(f,t)e^ndf  (9) 

Jo 

expresses  each  element  output  as  an  integral  in  2D  normalized  frequency  /  =  [/i ,  /2]  over  a  unit  square, 
one  period  of  the  integrand.  Because  that  integrand  is  periodic  with  the  region  of  integration  as  a  period, 
and  because  that  region  has  unit  area,  inverse  transform  Q  can  also  be  written  as  just  the  average  value  of 
the  integrand,  averaging  over  vector  /,  i.e., 

s°n(t)  =  avg{nf,t)e^}. 

This  lets  us  do  what  amounts  to  a  change  of  variable  in  the  integral  without  fussing  with  determinants  of 
Jacobians,  because  those  determinants  would  simply  deal  with  changes  to  the  1/area  factor  that  expresses 
an  average  as  an  integral.  Using  an  average  keeps  those  messy  factors  swept  under  the  mathematical  rug. 

Let  us  in  fact  use  change  of  variable 

k  =  /B+  /  =  kB  (10) 


to  write 

sn(t)  =  avg{s°(kB,  t)  ej27rkBn} 

where  of  course  the  average  is  taken  over  the  array  plane  only,  as  that  is  the  range  of  k  in  transformations 
([Toll.  We  can  then  go  just  a  little  further  and  use  an  ordinary  inverse  Fourier  transform  in  continuous  time  to 
write  5°(kB,  t)  =  f  S°(kB,  f)  e-j27rft  df  so  that  the  above  average  becomes 


s°J t)  =avg<{y  S°(kB,  f)  e-i27r(ft+kBn)  dfj. 

That  this  represents  an  LTSI  system  with  output  spatially  sampled  by  evaluation  at  x  =  ABn  as  discussed 
in  Section  [3T]  above  can  now  be  made  explicit  using  two  canceling  factors  of  A  to  rewrite  the  above  as 

S°(kB,  f)  ej27r(ft+Akx)  df 


s°n(  t)  =  avg 

k 


x=ABn 


(ii) 
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The  LTSI  portion  of  this  is  interpreted  in  detail  in  the  next  section. 


If  embedded  filtering  were  real  with  real  outputs,  then  conjugate  symmetry  S°(kB,f)  =  (S°(— kB,  —  f))*, 
which  would  put  components  averaged/integrated  over  into  conjugate  pairs,  would  be  guaranteed.  But  we 
see  in  Section  3.4  below  that  we  should  eliminate  the  system’s  vulnerability  to  image  beams  by  removing 
half  of  each  such  pair  with  filtering,  and  it’s  convenient  to  refer  that  filtering  into  the  element  system.  With 
the  balance  between  complex  pairs  thus  disturbed,  the  element  outputs  (t)  can  no  longer  be  assumed  real. 
They  must  be  complex. 


3.3  Physical  Interpretation  of  the  Element-Output  Averaging  Operation 

3.3.1  A  Complex  Plane  Wave 

To  interpret  element-output  representation  (fTT|)  correctly,  we  need  to  take  a  short  detour  into  the  mathe¬ 
matical  representation  of  electromagnetic  plane  waves.  For  mathematical  convenience  let  us  use  this  slightly 
unusual  parameterization  of  a  complex  plane  wave: 


4ej2"(ft+iMa'x)  . 


(12) 


That  this  is  a  propagating  wave  is  not  difficult  to  work  out.  If  we  let  9( t,x)  be  the  phase  of  the  complex 
exponential  so  that  9( t,  x)  =  27r(ft  +  ^  sgn(f)rt  -x),  we  can  predict  At  into  the  future  by  looking  in  the 
u  direction  by  Ax  so  that  9(t,  x  +  wAx)  =  9( t  +  At,  x)  if  we  set  f  At  =  -c  sgn(f)Ax  or,  equivalently,  by 
choosing  Ax  =  cAt  sgn(f).  Certainly  we  could  also  include  a  component  normal  to  u  in  Ax,  but  it  would 
only  be  lost  in  the  dot  product.  So  including  any  such  component  would  have  us  look  further  away  to  get 
the  same  information.  When  f  >  0  then,  Ax/ At  =  c  implies  that  the  wave  is  coming  from  the  u  direction  at 
speed  c.  Likewise,  Ax / At  =  — c  when  f  <  0,  and  the  wave  is  moving  in  the  u  direction  at  speed  c.  Evolving 
through  one  full  phase  cycle  requires  At  =  l/|f|  or  Ax  =  c/f,  so  the  wavelength  is  |Ax|  =  c/|f|. 

A  notationally  modified  version  of  complex  plane  wave  (fT2|)  is  more  useful  to  us.  Suppose  we  decom¬ 
pose  vector  y/|y|  u  into  a  sum  of  two  row  vectors,  an  array-plane  component  k  and  a  component  k  in  the 
boresight  direction.  Replacing  ^  1  f  u  in  (fT2])  with  y(k  +  k  )  yields  a  wave  in  the  form 


AejXk^x  ej27r(ft+Akx) 


(13) 


The  length  of  k  is  restricted  of  course,  because  ||k  +  kj_||  =  A||A(k  +  k±)||  =  A||^yit||  =  ^Ay.  The 
orthogonality  of  k  and  k  then  yields  Pythagorean  relationship 


|k||2+||kj2  =  llk  +  k. 


2 


from  which  it  follows  that  ||k||  <  -Ay  with  ||kx||  fully  determined  by  ||k||  and  with  k±  itself  extending 
distance  ||kj_||  along  either  the  boresight  direction  or  its  opposite. 
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3.3.2  Our  LTSI  System  Representation  Refers  to  Plane  Waves 


Choose  any  value  for  vector  k  with  ||k||  <  .  If  ||k||  =  it  follows  that  k±  =  0.  If  ||k||  <  — ^ 

however,  there  are  two  possible  k±  values,  and  they  are  negatives  of  each  other.  This  k  then  is  associated 
with  either  one  or  two  possible  incident  electromagnetic  plane  waves.  Such  an  incident  electromagnetic 
wave  has  x,  y,  and  z  electric-field  and  magnetic-field  components,  and  each  of  these  six  components  takes 
the  form  of  a  conjugate  pair  of  complex  plane  waves.  One  half  of  the  pair  takes  form  (|T3j)  with  f  >  0,  and, 
to  make  the  sum  of  the  two  halves  of  the  pair  real,  the  second  half  is  just  the  conjugate  of  the  first.  That 
conjugate  also  takes  form  (12 1  and  simply  has  different  A,  f,  k,  and  k  parameters.  In  particular,  because 
we  assume  the  first  half  of  the  conjugate  pair  had  f  >  0,  the  second  half  must  have  f  <  0. 


Consider  the  f  <  0  waves  of  these  field  components  in  more  detail.  There  are  6  such  waves  if  ||k||  =  yyy 
and  12  of  them  if  ||k||  <  ^y.  If  we  use  i  =  1, . . . ,  ?max  with  zmax  either  6  or  12  to  number  them,  we  can  say 
that  the  zlh  such  component  wave  takes  the  form 


Ci(k,f)ej2"(ft+ikx)  (14) 

with  one  value  of  f  and  one  value  of  k  common  to  all  such  component  waves  and  with  complex  amplitude 
Cj(k,  f)  corresponding  to  factor  A  el  if k  x  in  complex  wave  (fj~3|).  Among  our  zmax  component  waves,  only 
these  complex  amplitudes  differ,  even  having  different  units  according  as  the  component  represented  is  an 
electric  or  magnetic  field.  The  relationships  among  these  complex  amplitudes  are  determined  by  Maxwell’s 
equations  in  the  usual  ways,  but  none  of  that  concerns  us  here.  We  can  more  or  less  ignore  that  level  of 
detail  and  just  consider  these  basics,  which  are  true  no  matter  what  the  polarization  or  other  properties  of 
the  physical  wave  might  be. 

Our  LTSI  system,  being  linear,  responds  to  these  components  together  as  the  sum  of  its  responses  to  them 
separately^  Further,  an  LTSI  system’s  response  to  a  complex  exponential  in  t  and  x  is  that  same  complex 
exponential  scaled  by  a  frequency  response,  here  a  function  of  temporal  frequency  f,  spatial  frequency  k, 
and  of  course  the  k  and  component-identification  information  associated  with  i.  So  if  we  let  Gf(k.  f)  be 
that  frequency  response  for  the  zth  component,  our  LTSI  system’s  output  becomes 


‘■max  /  \ 

X]  Gi( k,  f)  Ci( k,  f)  ejMft+^kx)  .  (15) 


We  can  just  as  well  invoke  linearity  again  (and,  technically  speaking,  an  appropriate  form  of  bounded- 
input,  bounded-output  stability  as  well)  and  suppose  that  each  incident  field  component  actually  comprises 
an  integral  combination  of  corresponding  components  from  many  such  incident  waves  so  that  the  output 
becomes 

avgj  [  X^Gi(k,f)C'j(k,f)ej27r(ft+Akx)  df 

k  i= l 

*Of  course,  Maxwell's  equations  tell  us  that  the  imax  complex  amplitudes  associated  with  one  physical  wave  cannot  be  deter¬ 
mined  independently,  and  as  a  result  we  cannot  unambiguously  determine  imax  separate  response  characteristics.  But  that  does  not 
really  affect  our  argument  here.  There  would  certainly  be  no  harm  in  reducing  the  number  of  terms  to  account  for  this  lowering  of 
rank,  but  the  rest  of  the  argument  would  be  essentially  the  same  either  way. 
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Here  f  has  been  allowed  to  range  over  negative  as  well  as  positive  values  so  that  both  halves  of  conjugate 
pairs  of  waves  are  included.  Comparing  this  to  the  pre-sampling  average  in  element-output  representation 
CD.  we  see  that  for  k  with  ||k||  <  we  can  interpret  S°(kB,  f)  as  Y^t=\  ^(k,  f)  C',(k,  f).  Therefore,  for 
such  k  the  quantity  S°(kB,  f)  characterizes  the  LTSI  system’s  total  response  to  all  incident  complex  waves, 
of  whatever  polarization,  that  have  frequency  f  and  DOA  parameterized  by  k. 


3.3.3  The  Visible  Region  and  Choice  of  X 


For  a  narrowband  array  it  is  natural  to  choose  our  normalizing  distance  A  as  wavelength  c/|f|  so  that 
=  1,  so  that  our  k  for  plane  waves  is  just  the  array-plane  component  of  DOA  unit  vector  u.  This 
choice  results  in  easy-to-remember  bound  ||k||  <  1.  For  a  wideband  or  tunable  array,  however,  if  A  is  truly  a 
constant,  as  is  our  intent,  we  cannot  actually  have  A  equal  to  wavelength  c/|f|  for  more  than  one  frequency. 
Grating-lobe  considerations  discussed  below  in  Section  |5.1.3  make  it  convenient  in  this  case  to  choose 
A  =  c/|fj  for  |f|  at  the  top  end  of  the  band  of  interest.  At  that  frequency  then,  k  becomes  the  array-plane 
component  of  u  and  for  plane  waves  ranges  over  the  closed  unit  disk,  the  unit  circle  and  the  area  within  it. 
At  lower  frequencies  our  requirement  that  k  be  the  array-plane  component  of  u  means  k  ranges  over  a 
smaller  disk.  In  all  cases  the  range  of  k  corresponding  to  plane  waves  in  element  output  (fTT|)  is  called  the 
visible  region  of  the  k  plane  or  just  visible  space.  The  rest  of  the  plane  is  then  termed  invisible  space.  A 
specific  value  of  k  is  said  to  be  visible  or  invisible  according  to  the  region  in  which  it  lies. 


An  invisible  k  by  definition  does  not  correspond  to  any  plane  wave.  But  if  the  element  system  is  in  the 
far  field  of  all  sources,  all  electromagnetic  fields  take  the  form  of  plane  waves,  and  in  the  average  in  element 
output  ([IT])  these  arc  fully  accounted  for  by  the  portion  of  the  region  of  integration  that  intersects  with  visible 
space.  In  a  far-field  setting  the  portion  of  the  integration  that  intersects  invisible  space  must  not  contribute 
to  the  average,  so  for  invisible  k  complex  amplitude  S°(kB,  f)  =  0.  More  generally,  invisible  components 
of  S°(kB,  f)  represent  near-field,  evanescent  field  components.  In  this  report,  however,  far-field  conditions 
as  given  and  the  possibility  of  near-field  components  is  ignored. 


3.4  Assume  Analytic  Filtering  and  Take  the  Element  Output  at  the  Origin  as  a  Reference 

General  element-output  form  (fTT|)  reduces  when  n  =  0  to 


s°(t)  =  avg  j  J  S°(kB,  f)  df  |  (16) 

=  avg{5°(kB,t)}  (17) 

k 

where  an  ordinary  inverse  Fourier  transform  has  been  replaced  by  an  appropriate  function  of  t.  From  here 
on  let  us  assume  that  the  filtering  referred  into  the  elements  from  downstream  signal  processing  includes 
analytic  or  positive -pass  filtering.  This  lets  us  assume  that  the  integrand  in  (fT6])  is  effectively  zero  for  f  <  0. 
Such  an  assumption  is  most  convenient  because  it  resolves  the  direction  of  arrival/propagation  question 
unambiguously:  k  is  always  the  array-plane  component  of  which  points  in  the  DOA,  and  never 

its  negative.  Equation  ([17])  then  simply  expresses  the  output  of  the  element  at  the  origin  as  an  integral 
combination  of  components  originating  from  different  directions,  with  DOA  parameterized  by  k. 
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3.5  Relating  k  to  the  Wavenumber  Vector 

Readers  from  the  electromagnetics  community  might  naturally  wish  to  relate  their  familiar  notations  to 
what  is  used  here.  This  section  is  effectively  an  aside  to  address  that  and  can  be  skipped  without  harm. 

In  electromagnetics  a  complex  plane  wave  is  generally  written  J4ej(‘jt_kec'x*  with  radian  frequency 
uj  >  0.  The  “ec”  subscript  is  simply  to  distinguish  their  k  from  the  one  used  in  this  report.  The  negative- 
frequency  wave  of  a  conjugate  pair  is  seldom  written  explicitly  but  is  instead  generated  automatically  by  the 
conjugation  involved  in  taking  the  real  part  with  Re{z}  =  (z  +  z*) /2: 

Re{^4e^“t_kcc'x^}  =  |  e^“t_ktc'x^  -f  |V^li’t“kcc'x^  . 


The  two  complex  waves  are  expressed  in  terms  of  the  same  two  parameters  uj  and  kec  but  differ  by  a  sign 
on  the  entire  exponent.  By  contrast,  in  the  Section  3.3.2  discussion  above  the  two  complex  waves  of  a  pair 
presumably  look  alike  algebraically  but  simply  use  different  f  and  k  values.  Their  complex  amplitudes  also 
differ  from  those  used  here  by  a  factor  of  two,  but  this  is  really  irrelevant  since  we  never  refer  to  those 
amplitudes  alone  in  the  mathematics  but  only  to  products  representing  element  output  signals. 


More  important  is  how  direction  of  arrival  (DOA)  is  parameterized  in  the  two  notational  systems.  For 
simplicity  let’s  assume  f  >  0  and  compare  the  positive-frequency  waves.  Their  wave  above  matches  complex 
plane  wave  <[T2|)  when  u  =  27rf  and  kec  =  —  jy  jj  u ,  so  that  their  kec  points  in  the  direction  of  propagation  and 
has  length  equal  to  27t  over  wavelength.  For  /  >  0  our  k  is  the  array -plane  component  of  u ,  which 
makes  it  the  array -plane  component  of  —  ^  kec  ■  Our  DOA  parameter  k  can  be  obtained  from  their  radian 
wavenumber  vector  kec  by 


1.  multiplying  by  to  change  its  units  from  radians  per  unit  length  to  cycles  per  unit  length, 

2.  negating  it  to  change  its  direction  from  that  of  propagation  to  that  of  arrival, 

3.  multiplying  it  by  normalizing  constant  A  to  make  it  a  unit  vector  at  the  top  of  the  signal  band,  and 

4.  taking  the  array-plane  component. 


These  vectors  can  be  related  to  the  orientation  of  the  array  using  the  orthonormal  plotting  matrix  P  con¬ 
structed  above  in  Section  2.2.2  with  unit- vector  columns  in  azimuth,  elevation,  and  boresight  directions.  The 
antenna  community  typically  uses  direction  cosines  u,  v,  and  w  to  denote  components  of  DOA  unit  vector 
u  in  the  azimuth,  elevation,  and  boresight  directions  respectively.  Using  those  variables,  and  remembering 
that  k  is  a  row  vector  while  kec  is  typically  written  as  a  column  vector, 


■  U  ■ 

■  U  - 

u  =  P 

V 

k  —  —  22LP 

’  Kec  “  c/f 

V 

-W  - 

-W  - 

k  =  f{[uv  0]PT. 


For  planar  arrays  the  electromagnetics  community  typically  plots  “in  sine  space”  by  letting  horizontal  and 
vertical  axes  represent  u  and  v  respectively  and  using  using  either  a  color  scale  or  contours  to  represent 
DOA-dependent  function  values.  Because  u  is  a  unit  vector,  w  is  implied  by  u  and  v  and  need  not  be  shown. 
In  this  report  plots  assume  /  >  0,  in  respect  of  the  analytic-filtering  assumption  above  in  Section  3.4  and 
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horizontal  and  vertical  axes  represent  the  azimuth  and  elevation  components  of  k  respectively,  in  other  words 
the  first  two  components  of  kP  =  [it  v  0] ,  with  color  then  representing  function  values.  Therefore,  in 
the  usual  case  of  a  narrowband  array  with  A  set  to  c/|/|,  our  plots  are  the  same  as  those  of  the  antenna 
community. 


4,  MULTIDIMENSIONAL  DSP  BASICS 

We  need  very  little  from  multidimensional  DSP  actually,  just  a  very  few  basics. 

4.1  Three  Operators  and  a  Frequency-Referral  Identity 

The  element  output  indexed  by  n  is  written  everywhere  above  as  (t),  but  beginning  here  it  is  just  s°n. 
Our  processing  from  here  on  is  purely  spatial,  so  even  though  all  of  our  signals  remain  time  dependent,  the 
t  dependence  is  omitted  from  the  notation.  Just  s°,  with  no  subscript,  refers  to  the  entire  signal,  the  function 
of  n  with  all  its  many  spatial  samples,  and  that  concise  approach  is  used  for  other  discrete-space  signals  as 
well. 

Three  simple  discrete-space  operators  are  central  to  phased  arrays.  The  first  is  of  course  quite  standard. 
The  other  two  are  not,  but  they  are  most  convenient  when  analyzing  systems  with  phase-shift  steering. 


convolution 

s-kh  of  discrete-space 

quantities  s  and  h 


(s-kh)n  4  Ysmhn. 


phase  delay 
S  rv  /a  of  signal  s 
by  phase  gradient  fA 

steering 

/a  rv  h  of  impulse  response  h 
by  frequency  fA 


(sr,fA)n^sne-^n 
(fA^h)n  4  e/2^n  hn 


(18) 

(19) 

(20) 


Our  order-of-operations  convention  is  that  the  three  operations  above  are  done  after  multiplication  and  di¬ 
vision  but  before  addition  and  subtraction.  Unparenthesized  sequences  of  *  and  rv  operations  are  evalu¬ 
ated  left  to  right,  so  that  sr\fAkh  means  ( srxfA )  */i,  and  unparenthesized  sequences  of  *  and  ^  opera¬ 
tions  are  evaluated — this  is  not  standard  but  is  most  convenient — right  to  left,  so  that  h 1  *  fA  r~  lr  means 

/iM/a^  h2)- 

Recall  that  when  ordinary  filtering  after  a  frequency  conversion  is  to  be  replaced  by  filtering  before  the 
conversion,  the  filter’s  frequency  response  has  to  be  modified  appropriately.  We  can  derive  the  equivalent 
idea  for  spatial  processing  using  very  simple  steps: 


(s  rv  fA*h)n  =  ((S  rv  /a)  *  h)  n  (order  of  operations) 

=  rxfA)n—m  hm  (convolution  definition  (p~8]>) 

m 

=  sn-m  e-l27r^n-m)  hm  (phase-delay  definition  (Th])) 

m 
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—  (x^  *  pj27 r/Am  >  \  -j27r/Ai 

—  I  /  J  ^n—rri  c  '  ,jm  J  c 

rn 

=  e~y2nfAn 

rn 

=  {sk(fA^h)  r\fA)  n 


(split  the  complex  exponential) 

(steering  definition  (f20|)) 

(convolution  definition  (p~8]>) 
(phase-delay  definition  ([T9|)). 


This  frequency-referral  identity  takes  an  easy-to-remember  form  when  written  using  whole-signal  notation 
with  one  line  above  the  other  as 


S  rv  fA  k  h 

=  S  k  h )  fA.  (21) 


The  mnemonic  convenience  of  having  the  operands  in  the  same  order  on  the  page  on  both  sides  is  the  real 
reason  for  the  nonstandard  right-to-left  convention  adopted  for  sequences  containing  the  operation. 

4.2  Three  Simple  Fourier  Properties 

Four  simple  Fourier  properties  are  important  in  this  report.  Each  is  analogous  to  a  common  ID  Fourier 
property,  but  they  are  proved  here  for  the  2D  case  as  a  review  for  the  rusty  and  to  increase  the  comfort  level 
of  readers  new  to  multidimensional  transforms. 


Average  of  transform:  Given  Fourier  pair  sn  the  average  of  the  transform  is  just  the 

value  of  the  sample  at  the  origin: 


so=  f  S(f)df  =  avg{5(/)}, 
Ju  f 


This  is  just  inverse-transform  definition  ([8])  with  n  =  0  and  with  the  integral  re-expressed  as  an  average. 


Convolution  property:  Space-domain  convolution  corresponds  to  frequency-domain  multipli¬ 
cation,  so  that  xn<-kX(f)  and  ynk^Y(f)  imply  (xky)n<-kX(f)Y(f). 


Fourier-transform  definition  (|7])  and  convolution  definition  ([18])  applied  to  zn  =  ( xky)n  together  yield 


Z(f)  =  J2(xky)ne~^fn 

n 

=  e"j27r/n 


n 


m 
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Multiplying  the  summand  on  the  right  by  unity,  in  the  form  of  canceling  product  e  ]2~frn  e J 27Tf'n  of  complex 
exponentials,  and  then  rearranging, 

Z (/)  =  E  E  x™  yn-m(e-Wm  e^m)  e~^n 

n  m 

-  V  r  P -j2?r/m  y-  p-j27T f(n-m) 

—  /  J  c  v  yn—m  c 

m  n 

Let  n'  =  n  —  rn  i n  the  inner  sum.  Summing  over  all  n  and  summing  over  all  n'  are  equivalent,  as  the  same 
terms  get  summed,  so 

Z(f)  =  £  e~i2nfm  £  yn,  e-j2-/”'  =  X(/)  y(/). 

m  n' 


Frequency  shift:  Fourier  pair  hn  H(f)  implies  Fourier  pair  ( r-,  h )  n  <->  H(f-fA)- 


Fourier-transform  definition  ([7])  gives  H(f)  =  hn  e  J 2nfn.  Replace  f  with  /  —  /A to  obtain 

n 

=  J2(hne^)e-^n 


=  J2^h)ne~^fn 

n 


using  steering  definition  ([20])  at  the  end. 
fA^h. 


The  last  line  is  just  Fourier-transform  definition  ([7])  applied  to 


Periodicity:  Every  2D  discrete- space  Fourier  transform  H{f)  is  periodic  in  the  sense  that 
H (/)  =  H{f-  k)  for  all  k  £  g2. 


Just  replace  /  with  /—  k  in  Fourier-transform  definition  H(f)  =  hn  e  ]2llfn  to  see  that 


H{f-k)  =  YJh  »e_i 


-j27r  (f-k)n 


=  '$2  hn 

n 

=  H(f ) 


using  the  fact  that  e J27r^any  mte£er)  =  \  and  that  inner  product  kn  is  an  integer  because  the  elements  of  vectors 
k  and  n  are  integers. 
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5.  CREATING  AND  STEERING  BEAMS 
5.1  The  Single-Beam  Output 


5.1.1  Array  Steering,  Weighting,  and  Combining  Creates  a  Beam 


Assuming  for  the  moment  that  only  a  single  beam  is  required,  the  array  output  is  created  from  the 
element  outputs  using  a  row  two-vector  phase  gradient  fA  and  using  combining  weights  hn  that  here  take 
the  form  of  a  discrete-space  impulse  response.  We  represent  that  array  output  here  as  the  n  =  0  sample  of 
a  fictional  signal  sn,  fictional  because  no  other  sample  is  ever  realized.  Here  the  first  line  specifies  the  array 
output  exactly  as  it  is  ultimately  implemented  from  definitions  ( J_8 )  and  (p~9l),  and  the  other  lines  begin  the 
analysis  of  why  it  works: 


/a*  h) o 

(22) 

(/a^  h)  r*,fA)  o 

(/a^  h))0. 

(23) 

The  second  equality  follows  from  frequency -referral  identity  ([2Tj),  and  the  third  equality  just  recognizes  that, 
by  delay  definition  (19]),  a  phase  delay  leaves  the  zero  sample  unchanged. 

Next  we  take  Fourier  transforms.  Our  frequency-domain  notation  assumes  Fourier  pairs 

hn  H(f) 


and  the  frequency-shift,  convolution,  and  average-of-transform  properties  of  Section [4~2] above  then  respec¬ 
tively  yield 


{fA^h)n  ^  H(f-fA ), 

S°(f)H(f-fA), 

(s°*(fA*sh))  o  =  avg  {S°(f)H(f-fA)}. 


The  last  of  these  allows  us  to 
defined  in  transformation  (fTO]). 


rewrite  array  output  (|23j)  with  DOA  parameterized  either  by  /  or  by  k  as 
In  the  latter  case  we  also  need  analogous  relationships 


kA  =/aB+, 


/a  =  kAB. 


(24) 


The  two  forms  of  our  result  then  are 


so  =  avg{5(  f  )H(  f  -fA  )} 
/ 

=  avg{S(kB)tf((k-kA)B)}. 
k 


(25) 
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The  latter  is  just  zero-element  output  s{j(t)  from  ([17])  but  now  with  i?((k—  kA)B),  termed  the  array  factor, 
scaling  the  net  complex  amplitude  5(kB)  of  any  signal  components  arriving  from  the  DOA  specified  by  k. 
Unsteered  array  factor  H  (kB)  is  designed  to  have  a  narrow  beam  in  the  k  =  0  boresight  direction,  and  the 
beam  is  then  steered  by  setting  kA  to  the  value  of  DOA  parameter  k  where  the  beam  is  desired.  For  this 
reason  array-plane  three-vector  kA  might  reasonably  be  termed  a  steering  vector. 

5.1.2  Array-Factor  Periodicity 

Array  factors  are  periodic.  To  see  this,  begin  with  unsteered  array  factor  H (kB)  and  offset  DOA  pa¬ 
rameter  k  by  an  arbitrary  point  from  array-factor  periodicity  lattice  Z2B  f,  the  dual  of  the  element-position 
lattice.  For  an  arbitrary  kcZ2  then, 

H(( k  -  fcB+)B)  =  H (kB  -  feB+B) 

=  H (kB  -  k) 

=  H  (kB). 

The  second  and  third  equalities  follow  respectively  from  identity  B+B  =  I  and  the  periodicity  property  of 
the  Fourier  transform  in  Section  l4~2l  above. 

5. 1.3  Grating-Lobe  Prevention  Determines  Element  Spacing 

On  the  upper  left  and  upper  right  respectively  of  Fig.  [6]  are  an  example  element-position  lattice  and  the 
corresponding  array-periodicity  lattice,  respectively  comprising  points  in  the  normalized  position  plane  and 
in  the  beamspace  plane,  the  latter  a  function  of  DOA  parameter  k.  A  visual  representation  of  array-factor 
periodicity  lattice  Z2B+  like  that  shown  can  also  be  taken  as  a  sort  of  schematic  of  the  ideal  unsteered 
array  factor  itself.  This  is  because  unsteered  array  factor  II (kB)  is  most  commonly  designed  to  have  a 
narrow  beam  at  boresight  or  k  =  0,  a  point  of  array -periodicity  lattice  Z2B+.  The  periodicity  of  array  factor 
H (kB)  then  gives  it  a  beam  at  every  other  point  of  that  lattice  as  well,  so  that  all  lattice  points  schematically 
represent  unsteered  array-factor  beams. 

As  discussed  above  in  Section [373]  a  propagating  plane  wave  puts  k  in  the  visible  region  of  beamspace. 
There  it  is  just  the  array -plane  component  of  where  unit  vector  u  points  in  the  DOA.  When  frequency 

|f|  is  at  the  top  of  the  signal  band,  scale  factor  A|f|/c  becomes  unity  and  makes  k  just  the  array -plane 
component  of  the  DOA  unit  vector  u  itself.  When  we  draw  or  plot  the  beamspace  k  plane  then,  we  generally 
overlay  a  visible -region  globe  gridded  with  latitude  and  longitude  lines  representing  azimuth  and  elevation 
angles  of  this  DOA  unit  vector.  We  draw  it  with  unit  radius  to  represent  the  top  of  the  signal  band  and  simply 
understand  that  it  shrinks  proportionally  with  frequency  elsewhere.  When  unsteered  array  factor  H( kB)  is 
translated  or  slid  in  k  by  steering  vector  kA  to  create  steered  array  factor  II ((k  —  kA)B),  the  boresight 
beam  moves  to  kA  and  we  can  read  the  corresponding  azimuth  and  elevation  angles  from  the  visible-region 
globe’s  grid.  We  seldom  make  this  translation  explicit  in  the  graphics  and  instead  leave  it  to  the  imagination. 

In  the  beamspace  k  plane  on  the  upper  right  of  Fig.  [6]  however,  to  illustrate  a  point,  one  translated 
version  of  the  lattice  is  shown  very  with  the  schematic  boresight  beam  translated  to  —65°  azimuth 

and  —55°  elevation.  When  the  boresight  beam  moves  to  these  coordinates,  all  the  other  beams,  depicted 
prior  to  steering  by  lattice  points,  move  in  parallel.  In  the  example  one  such  steered  schematic  beam  appears 
just  to  the  upper  right  of  the  visible -region  globe. 
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Fig.  6  —  Top:  element-position  lattice  B$2  (left)  with  classic  nearest-neighbor  spacing  l/-\/3 
and  corresponding  array-factor-periodicity  lattice  Z2B+  (right)  with,  to  ensure  no  grating  lobes, 
a  nearest-neighbor  spacing  of  2,  the  latter  lattice  shown  both  untranslated  and  to  move 

the  boresight  point  to  an  example  beam  position  at  —65°  azimuth  and  —55°  elevation.  Bottom: 
closeups  of  Fig. [8j(left)  and  Fig.|7](right).  The  left  plots,  both  of  normalized  position  ^x,  are  on  the 
same  scale.  The  right  plots  show  the  k  plane  with  5°  grid  lines  and  with  the  lower  plot  zoomed  12  x 
relative  to  the  upper  plot. 
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Here  we  see  why  the  classic  choice  for  nearest-neighbor  spacing  in  array-periodicity  lattice  Z2B  +  is 
exactly  2.  In  array  design  increasing  element  density  B  1  B|~  1  //2/ A2  in  the  array  plane  increases  cost, 
so  generally  that  density  is  minimized.  Equivalently,  the  beamspace  density  B  1  Bj 1  2 A2  of  array-factor 
periodicity  lattice  Z2B+  is  maximized.  That  maximization  must,  however,  never  allow  array-factor  steering 
to  put  two  beams  in  the  visible -region  globe  simultaneously,  for  then  one  of  those  two  beams  would  be  an 
undesired  grating  lobe ,  potentially  causing  the  array  to  admit  contaminating  signals  from  some  undesired 
direction.  The  minimum  safe  spacing  for  these  periodic  beams  is  the  globe’s  diameter  of  2,  assuming  the 
boresight  beam  could  be  steered  until  its  edge  just  touched  the  globe  perimeter’s  inner  edge.  (Restrictions 
on  steering  in  some  systems  permit  nearest-neighbor  spacings  of  less  than  2.)  This  argument  is  for  the 
upper  end  of  the  signal  band,  but  the  globe  gets  smaller  for  lower  frequencies,  so  a  spacing  that’s  grating- 
lobe  safe  at  the  upper  end  is  grating-lobe  safe  at  lower  frequencies  as  well.  Given  such  a  lower  bound  on 
nearest-neighbor  distance,  the  beamspace  density  of  array-periodicity  lattice  Z2B  '  is  maximized  by  making 
it  a  hexagonal  lattice.  This  makes  element-position  lattice  Bzjc2  hexagonal  as  well,  and  the  spacing  of  2  in 
Z2B+  results  in  a  spacing  of  1  /yd  in  B^2,  giving  the  physical  element  centers  the  classic  spacincQof  A/\/3- 

5.1.4  Single-Beam  Implementation 


To  implement  array-output  computation  (|22|),  first  substitute  the  right  half  of  change  of  variable  ([24  > 
into  phase-delay  definition  (l9[)  to  obtain 


(s°r^fA)n 


p-j27rkABn 


and  then  use  convolution  definition  ( 18 1  to  see  that 


fs°rv 


/a*  fc)„  =  £«.e"I2"kaBm)  K-, 

m 


To  complete  the  implementation  of  array  output  ( 22 1,  just  sample  at  n  =  0  to  obtain 


«o  = 


£  h 


—  777.  I  5 


0  -j27TkABra^ 


o°  ^p-j27rkABra 

-m 


(26) 

(27) 


Computation  order  ([26])  is  typically  used  in  analog 
equally  reasonable  for  digital  arrays. 

5.2  Steering  Vectors  from  a  Steering  Lattice 


arrays,  and  either  computation  order,  (26)  or  (27), 


is 


The  single-beam  approach  above  can  of  course  be  used  several  times  in  parallel,  with  different  steering 
vectors  kA,  to  obtain  several  beams  simultaneously.  When  more  than  a  very  few  beams  are  to  be  realized 
digitally,  however,  an  FFT  approach  to  steering  is  likely  to  require  less  real-time  computation  in  the  imple¬ 
mentation.  FFT  steering’s  one  limitation  is  that  all  steering  vectors  kA  must  be  restricted  to  some  steering 

1  When  elements  are  required  to  be  on  a  square  lattice  in  spite  of  its  nonoptimality,  a  spacing  of  2  in  the  array-periodicity  lattice 
results  instead  in  the  equally  classic  A/2  physical  element  spacing  and  an  element  density  higher  than  in  the  A/%/3  hexagonal-lattice 

case  by  a  factor  of  2/%/3,  an  element  density  penalty  of  some  15.5%. 
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lattice  JE2Rste1erB+.  Because  Rsteer  must  be  a  2  x  2  resampling  matrix,  certain  sublattice  and  superlattice 
relationships  are  necessary. 


element-position  lattice 
dual  steering  lattice 


B%2  Z2B+ 

U  density  ratio  |  Rsteer  |  fl 

BRsteer?2  4 - 4  ^2ftsteer-E5  + 


array-factor  periodicity  lattice 
steering  lattice 


To  determine  the  effect  of  this  steering-lattice  restriction  on  array-output  computation  ([27]),  set  kA  = 
fcR-iB  1  for  some  k€^,2  and,  simply  to  put  the  result  into  the  most  familiar  form,  change  mton: 


so 

beam  A; 


^2(h_nS°n)  e-j2fffcR^B+Bn 

n 

J2(h-nS°n)  e-i2ri^"  . 


(28) 


The  second  equality  follows  of  course  from  identity  B+B  =  I. 

The  big -picture  diagram  of  beamspace  k  on  the  upper  right  of  Fig.  [6]  is  duplicated  in  Fig.  [7] in  a  much 
larger  size  with  several  additions,  most  notably  an  2rRjlL!clB  f  based  on  Rsteer  = 

[604  6°4] .  Each  point  of  that  is  enclosed  in  a  tiny  diamond-shaped  tile  for  increased  visibility. 

Of  far  more  significance  arc  the  large  tiles,  which  are  discussed  below  in  Section [574]  The  finer  features  of 
this  diagram  are  more  easily  visible  in  the  close  up  on  the  lower  right  in  Fig.  [6]  which  is  at  an  approximately 
6x  magnified  scale  on  the  page.  The  particularly  simple  form  of  resampling  matrix  Rsteer  used  here,  just  an 
identity  matrix  times  a  scale  factor,  makes  steering  lattice  S2R^e1erB+ just  a  scaled-down  version  of  array- 
factor  periodicity  lattice  Z2B+.  The  analysis  of  Section  5.5  below  implies  that  such  an  Rsteer  with  the  scale 


factor  a  power  of  two  leads  to  the  simplest  possible  structure  for  the  beamsteering  FFT. 


5.3  The  DFT  Input  Space 


Making  the  steering  lattice  a  superlattice  of  the  array-factor  periodicity  lattice  gives  identical  phases  to 
the  complex  exponentials  in  many  of  the  terms  of  multibeam  array-output  computation  (|28]).  Working  out 
the  details  here  allows  us  to  factor  out  those  common  factors  and  organize  beam-computation  inputs  into 
DFT-  or  FFT-style  bins.  The  key  is  a  particular  decomposition  of  element  index  n. 

As  we  work  this  out,  we  can  use  as  an  example  the  system  that  we  began  discussing  with  Figs.  [6]  and 
[7]  The  array-factor  periodicity  lattice  Z2B+  pictured  in  Fig.  [7]  and  on  the  right  in  Figs.  [6]  is  the  dual  of 
the  element-position  lattice  B^2  shown  in  Fig.  [8]  and  shown  about  5x  magnified  through  a  smaller  viewing 
window  on  the  lower  left  in  Fig.  [6j  The  2,527  specific  points  of  element-position  lattice  Bf?  that  are 
highlighted  in  Fig.  |s|  a  few  of  which  appear  in  Fig.  |6|  as  well,  arc  those  points  Bn  on  or  inside  a  circle 

of  radius  \J 69 1  /3  ~  15.18  centered  on  the  origin.  In  this  hypothetical  design  the  outputs  of  the  elements 
at  these  positions  are  the  ones  used  to  compute  the  array  output:  the  output  of  the  element  at  such  a 
position  Bn  is  multiplied  in  multibeam  computation  (|28|)  by  a  nonzero  weight  h-n,  making  these  the  active 
elements.  There  is  nothing  special  about  this  array  size  or  shape,  but  the  ideas  are  easier  to  discuss  in  the 
context  of  a  specific  example. 
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Fig.  7  —  Tile  boundaries  partition  beamspace  into  large  DFT  output  tiles,  identical  except  for 
translation,  such  that  each  point  of  ^2Rs4e2rB+  lies  in  exactly  one  tile  and  such  that 

each  tile  contains  exactly  one  point  of  array-factor  periodicity  lattice  Z12B+.  In  this  example  Rsteei-  = 
[6q46°4]  ,  so  the  tiles  each  contain  642  points  of  Z2R^e2rB+.  Here  a  simple  choice  of 

tile  shape,  based  on  the  basis-vector  rows  of  B  ,  lines  up  64  points  of  the  latter  along  each  tile  edge. 

divide  each  tile  into  a  64  x  64  grid  simply  to  make  the  structure  easier  to 

see,  both  here  and  in  the  closeup  on  the  lower  right  in  Fig.  [6] 
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Fig.  8  —  partition  the  normalized  array  plane  into  large  DFT  input  tiles,  identical 

except  for  translation,  such  that  each  point  of  element-position  lattice  Belies  in  exactly  one  tile  and 
such  that  each  tile  contains  exactly  one  point  of  BRstCci^2.  Positions  of  actual 

active  array  elements,  those  that  are  given  nonzero  weights  in  the  array-output  computation,  are 
highlighted  .  In  this  example  Rsteer  =  [6Q4  6°4] ,  so  the  tiles  must  each  contain  642  points  of  B^2 
Here  a  simple  choice  of  tile  shape,  based  on  the  columns  of  BRstecr,  lines  up  64  points 

of  the  latter  along  each  tile  edge.  A  closeup  appears  on  the  lower  left  in  Fig. [6] 


36 


Jeffrey  O.  Coleman 


5.3.1  Input  Tiling  and  Index  Decomposition 

We  begin  by  tiling  the  element-position  plane,  as  in  the  example  of  Figs.  [8]  into 


(defining  condition) 

identical  “input  tiles”  such  that  every  point  of  element-position  lattice  B  'f?  (or  the  whole  array 
plane  in  fact)  is  in  exactly  one  tile  and  such  that  each  of  the  tiles  contains  exactly  one  point  of 
dual  steering  lattice  BRsteer^2. 


Let  us  call  that  one  point  of  dual  steering  lattice  BRsteer^2the  tile’s  anchor.  Here  “identical”  means  any  tile 
can  be  translated — slid  with  no  change  in  orientation — to  create  any  of  the  others.  It  is  convenient  to  use  the 
tile  anchored  at  the  origin  as  a  sort  of  prototype  tile  in  that  way,  so  let  us  call  that  tile  the  (input  prototile) .  As 
long  as  the  defining  condition  above  is  satisfied,  it  makes  no  actual  difference  what  the  tile  shape  is,  but  it  is 
convenient  to  standardize  in  Sections  5.3.2|and|5.3.3  below  on  a  simple  way  to  use  basis  vectors  to  construct 
parallelogram-shaped  tiles  like  the  ones  in  Fig.  [8]  Regardless  of  tile  shape  the  defining  condition  for  tiles 
above  results  in  |Rsteer|  element  positions  in  each  tile,  simply  because  that  is  the  ratio  of  the  densities  of  the 
two  lattices. 


Given  a  satisfactory  tiling,  we  next  decompose  the  element  index  n  using  the  tiles  as  a  guide.  The 
strategy  comes  from  the  twin  facts  that  (1)  every  element  position  is  in  a  tile  and  (2)  every  tile  is  just  the 
(input  prototile)  translated  so  as  to  move  its  anchor,  the  origin,  to  a  desired  position.  It  follows  from  these 
facts  that  any  point  of  element  position  lattice  Bf2  can  be  uniquely  decomposed  as  a  sum  of  an  anchor 
position  from  dual  steering  lattice  BRsteer^2  and  a  point  of  element-position  lattice  B,I2  from  inside  the 
(input  prototile),  an  intra-tile  offset.  In  other  words,  for  any  n  £  2f2  there  is  exactly  one  combination  of  an 
m  G  If)2  and  an  [to]  E  %2  with  B[ro]  G  (input  prototile)  such  that 


Bn  =  BRsteerm  +  B[ro].  (29) 

Multiply  on  the  left  by  B+  and  use  identity  B+B  =  I  to  write  equivalent  decomposition 


n  =  Rsteer  m  +  M- 


(30) 


Here  the  outer  box  in  notation  M  signifies  that  its  possible  values  are  restricted  to  some  finite  set  determined 
by  a  tiling. 

5.3.2  The  Simplest  Case 

The  most  straightforward  system  designs  have  Rsteer  =  [  Tf*  Ayx]  *C)I  some  integer  expansion  factor  Nex. 
The  design  of  Figs.  [6]  [7]  and[8]is  such  a  system  with  Nex  =  64.  In  such  systems  steering  lattice  Z2R((CIL.IB 
is  just  array-factor  periodicity  lattice  Z2B+  shrunk  by  a  factor  of  Arex,  and  dual  steering  lattice  BRsteer?2 
is  just  element-position  lattice  B^2  magnified  by  Nex.  Each  tile  contains  exactly  |Rsteer|  =  A7^  element 
positions.  Using  notation 


n  i 
n2 


m  = 


~m  i ' 
777-2  . 


M  = 


'Mi ' 
M2. 


n  = 
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index-vector  decomposition  (|30|)  becomes  the  pair  of  equations 


ni  =  Nexmi  +  Mi, 
n>2  =  Nexm2  +  M-2- 


In  this  simple  special  case  a  simple  tiling  strategy  that  always  works,  the  one  used  in  Figs.  [6] through [8]  is  to 
choose  the  (input  prototile)  to  include  element  positions 

{Bn  such  that  n  =  [ )))  ]  with  0  <  m  <  Nex  and  0  <  n2  <  iVex}  .  (31) 

This  choice  results  in 


or 


Mi  =  ni  mod  Nex, 
\n\-2  =  n2  mod  Nex, 
mi  =  [ni/Nex\, 
m2  =  [n2/Nex\ 


M 

rn 


"  ni  mod  ./Vex 
.  n2  mod  A^x 

[ni/Nex\ 

[ Tl2 /  ^VexJ 


The  approach  just  outlined  is  actually  the  general  approach  discussed  next  applied  to  the  special  case  of  an 
Rsteer  of  the  form  [;'{f  A?J . 

5.3.3  Basis-Vector  Tiles 

Alternatively,  the  of  the  dual  steering  lattice  BRsteer/fi2  as  pictured  in  Fig.  [8]can  be  used  in 

a  slightly  more  mathematical  tile  construction  that  is  perfectly  general,  that  works  for  any  resampling  matrix 
Rsteer.  regardless  of  its  form.  An  element  position  is  naturally  expressed  as  p  =  B  n  for  some  n  e  ZjZ/2,  but 
we  know  how  to  express  that  position  in  any  basis  we  like.  In  Fig.  [8]  we  can  see  that  if  we  express  it  as  a 
linear  combination  of  the  of  the  dual  steering  lattice  and  write  it  as  p  =  BRsleer'«)  for  some 

column  two-vector  w  =  [  ] ,  point  p  is  in  the  (input  prototile)  if  both  0  <  w\  <  1  and  0  <  w2  <  1  hold. 

We  can  adopt  this  as  the  definition  of  basis-vector  tiles,  which  we  can  make  our  standard  scheme.  In  effect 
we’ve  constructed  our  (input  prototile)  as  a  parallelogram  with  corner  points  given  by  the  origin,  the  two 
basis-vector  columns  of  BRsteer.  and  the  sum  of  those  two  basis  vectors.  The  specific  choices  of  the  four 
“<”  and  “<”  tests  above  define  this  (input  prototile)  to  include  two  of  the  four  parallelogram  edges  and  one 
of  the  four  parallelogram  comers  while  excluding  the  others,  choices  that  satisfy  the  defining  condition  for 
tiles  in  S ection  15 . 3 . 1 1  above . 

To  find  w  such  that  p  =  BRslee[  w,  left  multiply  the  latter  by  RstelerB+  and  use  identity  B+B  =  I 
to  show  that  w  =  R0C(.rB  p,  then  substitute  the  original  form  p  =  Bn  and  use  identity  B+B  =  I 
again  to  obtain  w  =  R^rn.  Once  we  adopt  Matlab’s  convention  that  the  floor  operation  that  rounds 
towards  —  oo,  denoted  mathematically  by  enclosure  in  |  J  brackets,  should  apply  to  vectors  elementwise, 
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conditions  0  <  w\  <  1  and  0  <  wo  <  1  can  be  jointly  expressed  as  just  w\  =  0.  Our  (input  prototile) 
therefore  contains  element  positions  Bn  for  which  RZeer^J  =  0. 

This  actually  gives  us  a  simple,  general  way  to  calculate  m  and  [n],  because  the  B|n]  term  in  lattice- 
point  decomposition  (|29|)  is,  by  its  definition  as  an  intra-tile  offset,  itself  a  position  in  the  (input  prototile) 
and  therefore  has,  using  our  new  (input  prototile)  definition,  |_Rsteer20j  =  0-  This  makes  it  trivial  to  solve 
equivalent  index- vector  decomposition  (|3Q|)  for  m:  just  left  multiply  by  and  observe  that,  since  adding 
an  integer  can  always  be  pulled  out  of  a  floor  operation  so  that  [n  +  x\  =  n  +  \_x\ , 

LRsteer«J  =  [rn  +  RSteer[n]J  =  m  +  [R^rlSlJ  =  m- 
In  principle  then,  obtaining  index- vector  decomposition  (|30j)  then  amounts  to  just  computing 

m  =  LRs"teernJ , 

M  =  n  -  Rsteer  7™ 


(32) 

(33) 


5.3.4  Making  Numerical  Computation  Robust 


As  written,  anchor-point  index  computation  (32)  is  vulnerable  to  numerical  errors.  For  example,  in 


Matlab  on  the  author’s  computer,  using  Rsteer  =  [j  5]  and  n  =  [_19]  and  computing  Rsle'ern  naively 
using  inv  (R)  *n  or  thoughtfully  using  R\n  yields,  respectively, 


Rsteer™  = 


-7 

-1-2 


-51 


or 


R  1  77  — 
-“'steer 


-7 

-1-2 


-52 


when  actually  Rsle’ern  =  [  _( ]  would  be  correct.  These  numerical  emors  in  the  51st  and  52nd  bits  to  the  right 
of  the  binary  point  cause  Matlab’s  floor  (  )  function  to  yield  LRjjeernJ  =  [l2]  instead  of  the  correct 
result  LRsteer™J  =  [ll]- 


We  can  rework  computation  (|32|)  to  remove  this  numerical  vulnerability.  Multiply  the  argument  of  the 
floor  operation  by  both  |Rsteer|  and  its  inverse  and  write 


m  = 


|RSi 


IR: 


■steer 


Rsteer)  ^ 


We  can  round  the  elements  of  the  product  in  the  parentheses  to  integers  to  remove  the  numerical  noise  of 
matrix  inversion  if  we  can  argue  that  this  product  must  actually  be  an  integer  matrix.  For  any  invertible 
matrix  R  in  fact,  |R|  R_1  =  ±  det(R)  R_1  =  ±  adj(R),  using  in  the  last  step  that  R_1  =  det|R^  adj(R). 
When  R  is  an  integer  matrix,  its  adjugate  matrix  adj(R)  is  an  integer  matrix  as  well,  as  therefore  are 
±adj(R)  and  |R|  R-1.  Now 


rn  = 


IR, 


■steer  | 


rOUnd(|Rsteer|Rsteer)n 


(34) 
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For  any  integer  m  and  positive  integer  N,  one  variety  of  integer  division  is  defined  by  m  /  N  A  [m/N\ 
and  can  be  computed  in  Matlab  (even  for  vector  m )  as  double  (idivide  (int64  (m)  ,N,  '  floor'  )  ) . 
It’s  clumsy  but  absolutely  safe.  The  quantity  in  large  parentheses  in  computation  ([34)  is  an  integer  vector, 
so  using  such  an  integer  division  of  every  element  of  that  parenthesized  vector  by  |Rsteer|  allows  us  to  write 


m  = 


round!  |R, 


steer 


R 


steer 


(35) 


and  be  fully  confident  that  there  are  no  numerical  issues.  Indeed,  test  computation  yields  |  Rstc11.rnJ  =  [  _1] 
for  the  example  above. 

An  alternate  approach,  based  on  solving  index-vector  decomposition  (f30|)  for  m  after  first  computing 
|n]  safely  using  Matlab’s  mod  (  )  function  rather  than  integer  division,  is  presented  in  Ref.  33  Integer 
matrix  ro  u  nd(  |  Rsleeij  R(teer)  is  also  at  the  heart  of  this  alternate  approach.  Regardless  of  which  approach  is 
taken,  this  integer  matrix  is  just  if  det(Rsteer)  >  0  or  [~rrf1  _™]  if  det(Rsteer)  <  0,  where 

tj  _  r  rn  ri2 1 

-TV-steer  —  [  r2 1  r22  J  ’ 

5.3.5  Folding  Weighted  Inputs  Creates  Input  Bins 


Since  decomposition  (f30|)  of  n  is  always  available,  in  multibeam  computation  (28 1  we  can  substitute  for 
n  and  sum  over  m  and  WV\  instead: 


■so 

beamfc 


—  ^  ^  G  J  ^Tr&R-steerlZll 

m  G  % 2  such  that 
B|n]  G  (input  prototile) 


(36) 


A 


E 


[h-nS°n] 


n= Rsteer  m+@' 


(37) 


A  complex-exponential  factor  e-j27rfcR.stcerRsteerm  _  e-j27rfcm  _  ^  |ias  disappeared.  Engineering  this  disap¬ 
pearance  is  in  fact  the  whole  point  of  decomposition  using  input  tiles,  because  that  decomposition  is  what 
makes  the  only  remaining  complex  exponential  independent  of  m  so  it  can  be  factored  out  of  the  sum  in 
([37])  to  leave  the  preliminary  folding  (sometimes  termed  data  turning )  of  weighted  element  outputs  there  in¬ 
dependent  of  beam  index  k.  No  matter  how  few  or  how  many  beams  we  compute,  we  compute  (|37|)  exactly 
I  Rsteer  |  times,  once  for  each  possible  value  of  \n\,  in  effect  obtaining  what  in  the  next  section  are  the  |Rsteer| 
distinct  input  bins  of  a  generalized  DFT. 


What’s  less  obvious  but  no  less  important  is  that  folding  sum  ([37]),  a  sum  over  all  tiles  that  is  calculated 
for  each  position  within  a  tile,  in  practical  designs  has  very  few  nonzero  terms,  often  just  one  or  even  none. 
This  is  the  case  for  our  example  design,  in  fact,  as  can  be  seen  in  Fig.  [8j  In  that  example,  for  each  intra¬ 
tile  offset  BgH  there  is  at  most  one  m  for  which  the  corresponding  point  in  Fig.  jsjhas  a  nonzero  weight 
h-(nst^m+m)  in  folding  sum  ([37). 


5.4  The  DFT  Output  Space  and  the  Beamsteering  DFT 


Beamspace,  as  pictured  in  the  example  of  Fig.  [7] and  its  closeup  in  the  lower  right  of  Fig.  [6j  can  be  tiled 
with  identical  output  tiles  (which  need  not  resemble  the  input  tiles)  such  that  every  tile  contains  exactly  one 
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point  of  array-factor  periodicity  lattice  Z2B+  and  every  point  of  steering  lattice  Z2R,st(!erB+  is  in  exactly 
one  tile. 

Suppose  an  allowable  output  tiling  is  in  place.  The  tile  containing  an  arbitrary  beam  position  fcR^e1erB+ 
also  contains  exactly  one  point  jB+of  the  array-factor  periodicity  lattice  (with  j  not  to  be  confused  with 
j  =  from  which  integer  row  two- vector  E  can  in  principle  be  solved  for  (practical  approaches  to 

follow)  to  obtain  unique  decomposition 

fcR-1erB+=  jB+  +  |Rste1erB+, 

expressing  the  beam  position  as  a  point  jB+  in  the  array-factor  periodicity  lattice  plus  an  intra-tile  offset 
!R(,'e|.B+  from  the  steering  lattice.  Multiply  through  on  the  right  by  BRsteer  and  use  B+B  =  I  to  obtain 
the  equivalent 

k  =  j  Rsteer  +  El-  (38) 


In  the  common  Rsteer  =  jy  ]  sPecial  case — it  applies  with  integer  expansion  factor  Nex  =  64  in  the 
example  of  Figs.[6]through[& — the  (output  prototile)  anchored  at  the  origin  includes  beam  positions 

|A:R((e1erB+such  that  k  =  [fci  with  0  <  k\  <Nex  and  0  <  k2  <-/Vex}  .  (39) 

More  generally  we  can  always  construct  basis-vector  tiles  by  assigning  to  the  (output  prototile)  the  beam 
positions  £:RSJL,IB+  for  which  /cR-OJa-J  =  0.  Then  index  decomposition  ([38])  amounts  to  computing 


j  L^-^steerJ  > 

\B  =  k  J  Rsteer- 


To  realize  anchor-point  index  computation  (|40])  robustly,  use 

3  —  (k  r°und(  Rsteer  Rsteer)  j  R-stccr 


(40) 


(41) 


An  alternate  approach  HI  solves  index- vector  decomposition  ([38])  for  j  after  first  computing  E  safely  using 
Matlab’s  mod  (  )  function  rather  than  integer  division.  Integer  matrix  round) |Rsteer|R(teer)  can  be  directly 


constructed  as  outlined  above,  at  the  end  of  Section [53)4]  if  desired. 


Apply  the  output  tiling  to  beam  computation  by  substituting  unique  decomposition  (|38|)  into  multibeam 


computation  (36): 


so  y 

beam  j  Rsteer + E  ^  [ H 


(42) 


=  x0e_j27r®Rstt“111 . 

[n]G  (input  prototile) 


(43) 
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Recognizing  e_j2wj'RsteerR'st“r0  =  e-j27rj@  —  \  |ias  simplified  the  right  side  of  (|43|),  making  it  in  fact  inde¬ 
pendent  of  j .  That  independence  reflects  that  the  phase  shifts  required  to  steer  unsteered  boresight  beam  to 
a  given  position  (jRsteer+  H)R^e1erB+  don’t  depend  on  j,  because  the  array  factor’s  periodicity  makes  any 
steering  increment  of  jB  '  irrelevant. 


The  sum  on  the  right  in  (43 1  is  a  generalized  2D  discrete  Fourier  transform  (DFT)  (2),  in  recognition 
of  which  the  traditional  DFT  output  name  Xg  is  used  alongside  a  functional  designation  as  such-and-such 
beam  of  sq- 


5.5  FFT  Realization  of  the  DFT  in  the  Simple  Case 

When  Rsteer  =  ['Tr  7vex]  f°r  some  integer  expansion  factor  Xex,  the  beamsteering  DFT  can  be  computed 
using  an  ordinary  2D  FFT.  This  is  developed  here,  but  this  simple  special  case  must  not  be  confused  with 
the  general  result.  More  generally  a  factorable  Rsteer  leads  to  the  generalized  FFT  structure  developed  later 
in  Section  1531 


In  the  simplest  case,  Rsteer  =  [Aqx  n  ]  f°r  some  integer  expansion  factor  Xex.  If  the  input  and  output 
prototiles  are  basis-vector  tiles  as  per  special-case  relationships  (|3T|)  and  ((39]),  we  can  use  expansions 


M  = 


E  =  [ii  i2], 


iRi®  =  (ISiiai  +  E2®2)/xex  (44) 


to  rewrite  DFT  ( 43 1  as 


tVgx  —  1 

0!,  @2  =  0 


n  p-j2ff(EiSi+g|2S2 

@i  c 

@2j 


JVex- 1  /tVex-l 

EE* 

@2=0  \@1=0 

Wx-l  /  tVex  1 

E  E* 

@1=0  \  @2=0 


mu 

@2 


@i 

@2 


3-j27rHi@i/tVex  1  p-j27r|J]2@2/Arc> 


x-j27r|i]2@2/Arex  I  p— j27r|fc]l@l/tVej 


(45) 

(46) 

(47) 


The  form  of  ([45])  is  that  of  an  ordinary  Nex  x  Nex  point  DFT  in  2D,  which  is  conventionally  computed  as  per 
either  of  equivalent  forms  ([46])  and  (|47]).  Figure  [9]  illustrates  special-case  DFT  (|46])  in  block-diagram  form. 
Each  vertical  block  on  the  left  realizes  the  inner  sum,  a  ID  DFT,  for  a  particular  value  of  [rah  to  transform 
the  dependency  on  position  index  [rrh  into  dependency  on  beamspace  index  \k\i-  Each  horizontal  block  in 
Fig.[9]then  realizes  the  outer  sum  for  one  value  of  Ei,  transforming  dependency  on  position  index  \n\2  into 
a  dependency  on  beamspace  index  E2-  Realizing  equivalent  special-case  DFT  ([47])  instead  would  transform 
first  on  \n\2  and  then  on  @1.  It  makes  no  difference  which  set  of  transforms  comes  first.  If  Nex  is  chosen 
to  be  highly  composite  in  the  usual  ways,  each  ID  DFT  block  in  this  structure  can  be  efficiently  realized 
with  an  FFT  algorithm,  and  it’s  both  appropriate  and  completely  conventional  to  then  consider  the  overall 
2D  realization  a  2D  FFT. 


In  many  radar  applications  only  a  cluster  of  beams  is  required,  and  many  of  the  beam  outputs  X\ 


[Ei  E2] 


are  unneeded.  This  typically  allows  some  output  FFTs  to  simply  be  removed,  and  of  course  the  input  FFTs 
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inputs 

•T\n\ 


Fig.  9  —  Example  simplest-case  FFT  digital  steering 
structure  for  Arcx  =  4.  Weighted  and  folded  element 
outputs  xg  are  here  processed  in  turn  by  a  bank  of  ID 
on  [n]i  and  a  bank  of  ID  FFTs  on  \n\2  to  create 
beam  outputs  Xgj .  The  components  of  index  vectors 
El  =  [  and  [0  =  [@1  i2]  range  over  0, . . . ,  Nex-1. 


then  need  not  provide  those  outputs  that  would  have  driven  the  removed  FFTs.  Ignoring  these  potential 
simplifications  in  favor  of  the  simplicity  of  using  FPGA  or  VLSI  “cores”  as-is  is  apt  to  be  expensive  both 
in  chip  area  and  power  consumption  and  may  in  borderline  cases  even  negate  the  advantage  of  the  FFT 
approach  over  independent  brute-force  computation  of  every  beam  of  the  cluster. 

5.6  FFT  Realization  of  the  DFT  in  the  General  Case 


Typical  reasons  to  consider  going  beyond  the  2D  FFT  just  discussed  are  outlined  next.  A  suitable  gen¬ 
eral  FFT  factoring  is  then  developed,  beginning  with  appropriate  tiling-based  element-  and  beam-position 
addressing  schemes. 


5.6.1  An  Extra  Small-Determinant  Factor  in  Rsteer  Can  be  Useful 


As  per  the  discussion  of  sublattice  density  at  the  end  of  Section  2.3.1  above,  there  are  |Rsteer|  beam 
positions  in  steering  lattice  Z2Rjjclc|.B^  per  period  of  array-factor  periodicity  lattice  Z2B+.  The  area  of  that 
period  is  generally  fixed  by  the  grating-lobe  considerations  discussed  above  in  Section  5.1.3  so  typically 
beam  density  is  actually  set  by  |Rsteer|  alone.  The  resampling  matrix  Rsteer  =  7vex]  of  the  last  section 
has  |  Rsteer  |  =  Ar2x,  so  when  Nex  is  limited  to  powers  of  two  changing  iVex  cannot  effect  changes  in  beam 
density  by  less  than  a  factor  of  four.  Design  options  for  |  Rsteer  |  then  start  at  some  initial  power  of  four,  the 
lowest  to  be  considered,  and  go  up  by  factors  of  four. 


4initialx  one  of 


40  41  44 

43  . 

= 

1  4  16 

32. 

Many  more  possibilities  result  from  allowing  a  little  more  flexibility  and  permitting  Nex  to  have  one 
factor  of  three  as  well  as  arbitrarily  many  factors  of  two,  as  one  factor  of  9  is  now  allowed  in  |Rsteer|. 


r 

4° 

41 

42 

43  44 . 

4initialx  one  of  \ 

x  1  = 

1 

4 

16 

64  256 . 

[  x  9  = 

9 

36 

144 

576 .. . 

In  effect  this  is  a  factor  of  9  or,  say,  16  times  an  arbitrary  power  of  four,  where  the  factor  of  9  is  the  new 
possibility  introduced  here.  A  density  change  of  16/9  or  9/16,  relative  to  densities  available  without  the 
factor  of  9,  becomes  an  option. 
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One  attractive  feature  of  the  Rsteer  =  [  Tf  tv  ]  f°rm  is  that  the  steering  lattice  has  the  same  geometry  as 
the  array- factor  periodicity  lattice.  Only  the  scale  is  changed.  A  hexagonal  array- factor  periodicity  lattice 
yields  a  hexagonal  steering  lattice,  and  a  square  array-factor  periodicity  lattice  yields  a  square  steering 
lattice.  If  we  allow  a  rotation  as  well  as  a  scaling  in  obtaining  the  steering  lattice  from  the  array- factor 
periodicity  lattice,  we  can  obtain  even  more  possible  | Rsteer |  values  while  otherwise  preserving  this  same- 
geometry  feature.  If  the  array-factor  periodicity  lattice  is  hexagonal,  for  example,  the  matrix  product 


R 


steer 


O) 

X 

0 

'1  -1' 

X 

0 

.1  2. 

(48) 


makes  the  steering  lattice  hexagonal  as  well,  but  scaled  and  rotated  30°.  Scaling  is  by 


R-steer  — 


O 

X 

£ 

1  -11 

X 

0 

.1  2J 

=  3  NL. 


This  optional  factor  of  three  increases  the  spectrum  of  |Rsteer|  possibilities. 
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In  effect  this  is  a  factor  of  9,  12,  or  16  times  an  arbitrary  power  of  four.  These  are  fairly  fine  steps. 

Instead  of  obtaining  a  factor  of  three  in  |Rsteer|  by  including  [  J  ”'],  we  could  obtain  a  factor  of  7  or 
13  or  19  by  including  [  2  ~ 4  ]  or  [ _4]  or  [  :J  “  l  ]  respectively,  with  different  resulting  steering-lattice 
rotations  in  each  case,  though  computational  efficiency  suffers  some  with  larger  determinant  factors,  even 
with  the  relatively  efficient  approach  of  the  next  section. 


If  the  array-factor  periodicity  lattice  were  square,  a  factor  of  2,  5,  13,  or  17  could  be  included  in  |Rsteer| 
by  respectively  including  in  Rsteer  =  [A(f  a?J  an  extra  factor  of  [  J  ~j],  [2  [jj  _3]’or  [i  "4]-  In 
each  case  the  steering  lattice  would  remain  square  but  would  be  rotated. 


A  systematic  approach  to  constructing  these  geometry-preserving  resampling-matrix  factors  is  described 
in  the  background  material  on  sublattices  in  Ref.  [6]  Other  possible  reasons  to  include  additional  matrix 
factors  in  Rsteer  are  discussed  in  Ref.  [3]  most  importantly  obtaining  a  nearly  square  steering  lattice  from  a 
hexagonal  array-factor  periodicity  lattice  and  vice  versa.  In  general  incorporating  a  new  small -determinant 
Rsteer  factor  should  be  viewed  as  an  available  tool  for  more  flexibly  customizing  the  steering  lattice. 


5.6.2  A11  Intermediate  Lattice 

Suppose  Rsteer  =  RoutRin>  for  example  with  one  of  Rout  and  Rm  of  the  form  [^ex  and  the  other 
a  resampling  matrix  with  a  small  determinant  as  discussed  in  the  last  section.  This  naturally  introduces  an 
intermediate  lattice  ^2R”ujB+  between  the  array-factor  periodicity  and  steering  lattices,  so  that  the  three 
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lattices  form  a  sublattice  chain,  as  do  their  duals. 


element-position  lattice 

B^2 

^2B  + 

array-factor  periodicity  lattice 

U  density  ratio  |  Rol,t  |  H 

dual  intermediate  lattice 

BRout|^ 

<— *•  22r0-u|b+ 

intermediate  lattice 

U  density  ratio  |  Rqn  |  H 

BRoutRin?2 

|| 

S2Rr1R-u{B+ 

|| 

dual  steering  lattice 

BRsteer^2 

< — >  ^  RsteerB 

steering  lattice 

5.6.3  Element  Addressing  Using  Two  Position-Space  Tilings 

We  can  use  the  lattices  on  the  left  in  relationships  (|49|)  to  develop  an  alternate  addressing  scheme  for 
element  positions  that  extends  the  position-space  tiling  ideas  developed  above,  in  Section [53] 

First,  the  basis  vectors  of  dual  intermediate  lattice  BR0Ut|2  can  be  used  to  construct  small  position- 
space  tiles,  as  pictured  in  the  example  of  Fig.  [TO}  to  tile  the  plane  such  that  each  tile  contains  one  point,  its 
anchor  point,  from  that  dual  intermediate  lattice  and  such  that  each  point  of  element-position  lattice  Rf? 
is  in  exactly  one  tile.  This  leads  to  a  unique  decomposition,  in  either  of  the  two  equivalent  forms,  of  an 
arbitrary  point  from  that  element-position  lattice  as  the  position  of  an  anchor  point  plus  an  intra-tile  offset, 

Bn  =  BRoutn'-f  B[n],  (50) 

n=  Routri/  +  m  (51) 


where,  if  we  use  the  optional  but  convenient  basis-vector-tile  approach  developed  above  in  Section  5.3.3 


n'=  (52) 

M=  n  —  Romn'. 


The  decomposition  can  be  determined  graphically,  using  a  diagram  like  that  of  Fig.  [TOj  or  the  approach 
developed  above  in  Section  5.3.4  can  be  used  to  compute  small-tile  anchor-point  index  ([52])  robustly.  Ei¬ 
ther  way,  the  values  that  can  be  taken  by  [n]  are  the  |R0Ut|  specific  integer  column  two- vectors  for  which 
[RouM  =  0. 


Second,  the  basis  vectors  of  dual  steering  lattice  BRsteer^2  can  be  used  in  a  similar  way  to  construct 
large  position-space  tiles,  as  pictured  in  the  example  of  Fig.  [TO}  to  tile  the  plane  such  that  each  tile  contains 
one  point,  its  anchor  point ,  from  that  dual  steering  lattice  and  such  that  each  point  of  dual  intermediate 
lattice  BRout|2  is  in  exactly  one  tile.  This  leads  to  a  unique  decomposition  of  an  arbitrary  point — point 
BRout'n/  from  decomposition  (|50[)  above  is  especially  interesting — from  that  dual  intermediate  lattice  as 
the  position  of  an  anchor  point  plus  an  intra-tile  offset  in  either  of  the  two  equivalent  forms 


BRoutn  —  BR0utRin^  T  BR.nntirn.|, 
n!  =  R;n£  +  [m] 


(53) 

(54) 
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Fig.  10  —  To  address  element  positions  with  three  digit  vectors,  use  the  Section|5.3.3|basis-vector 
tiling  twice  to  create,  in  the  example  here  with  Ri„  =  [  J  ]  and  Rout  =  1<y"j 

containing 

an  anchor  and  points 

point  from  exactly  from 

small  tiles  ?  ^  dual  intermediate  lattice  BRout^2  ^  'v  |  RolIt  |  ^  'v  element-position  lattice  B^2  and 
large  tiles  ^  dual  steering  lattice  BRoutRin^2  ^  jRm|  ^  dual  intermediate  lattice  BR0„t|2 
Each  element  position  can  then  be  uniquely  expressed  as 
a  large-tile  anchor  point  BRoutRin£ 

+  an  offset  BRout|m]  to  one  of  that  large  tile’s  small-tile  anchor  points 

+  an  offset  B[n]  to  one  of  that  small  tile’s  element  positions. 

The  (example)  has  l  =  [  _ }  ] ,  gn]  =  [  ?  ],  and  gg  =  [  g1  ].  Element  outputs  are  summed  into  the  FFT 
input  bins  in  the  |highTighted  region  located  by  setting  i  =  0  to  address  bins  with  gn]  and  gg  only. 
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where,  if  we  again  choose  to  use  the  basis-vector-tile  approach  of  Section  5.3.3  above 


t=  LRrVj, 

20  =  n' —  Rin-^. 


(55) 


The  decomposition  can  be  determined  graphically,  using  a  diagram  like  that  of  Fig.  [lOj  or  the  approach  of 
Section  5.3.4|can  be  used  to  compute  large-tile  anchor-point  index  (|55|)  robustly.  Either  way,  the  values  that 


can  be  taken  by  \m\  are  the  Rm  specific  integer  column  two-vectors  for  which  |  Rin  \fff  j  =  0. 


Decompositions  (50 )  and  ([5T])  can  be  combined  with  decompositions  (|53])  and  (54 )  to  yield  the  decom¬ 
positions  of  the  element  position  and  the  element -position  index  vector  given  by 


Bn  =  BRoutRinf  +  BKoutm  +  B[n],  (56) 

n  —  ROL,tRm£  T  Rout®  T  [ xf\-  (5 V) 


In  Fig.  10  the  three  thin  line  segments  leading  from  the  origin  to  the  one  (example)  element  position  cor¬ 
respond  to  the  three  terms  on  the  right  in  element-position  decomposition  ([56]).  (The  fourth  segment  is 
discussed  below.) 


An  element  position  is  indexed  or  addressed  by  column  vector  n  or,  using  address  decomposition  ([57]), 
the  three  column  vectors  i,  pm],  and  [70,  so  in  effect  (|57[)  relates  two  element-position  addressing  schemes. 
It  also  generalizes  on  mixed-radix  representation  of  an  integer.  Here  l,  [ml.  and  \n\  are  digit  vectors  from 
if2  and  can  respectively  take  an  infinite  number  of  values,  exactly  R,m  values,  and  exactly  |R0Ut|  values. 
Resampling  matrices  Rout  and  Rm  are  radix  matrices.  The  most  significant  vector  (MSV)  is  f,  and  the 
least  significant  vector  (LSV)  is  [70.  When  the  vectors  and  matrices  involved  are  replaced  with  integers, 
the  procedure  here  for  obtaining  the  digit  vectors  in  fact  degenerates  to  the  usual  procedure  for  obtaining  a 
mixed-radix  representation  of  an  integer. 

The  key  ideas  of  this  position-space  tiling  are  summarized  for  review  convenience  in  the  Fig.[l0]caption. 


5.6.4  Beam  Addressing  Using  Two  Beamspace  Tilings 

In  a  similar  way,  we  can  use  the  lattices  on  the  right  in  the  relationships  of  (|49[)  to  develop  an  alternate 
addressing  scheme  for  beam  positions  that  extends  the  beamspace  tiling  ideas  from  Section[54] above. 

First,  the  basis  vectors  of  intermediate  lattice  ^2R~u^B+can  be  used  to  construct  small  beamspace  tiles, 
as  pictured  in  the  example  of  Fig.  [TTJ  to  tile  the  plane  such  that  each  tile  contains  one  point,  its  anchor  point, 
from  that  intermediate  lattice  and  such  that  each  point  of  steering  lattice  Z2Rr1R“U(B+is  in  exactly  one 
tile.  This  leads  to  a  unique  decomposition  of  an  arbitrary  point  from  that  steering  lattice  as  the  position  of 
an  anchor  point  plus  an  intra-tile  offset  in  either  of  the  two  equivalent  forms 


^'Rou'tB+=  k'K  ~^B+  +  iR-1R-1|B+, 
k  =  fc/Rjn  +  M 


(58) 

(59) 
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Fig.  11  —  To  address  beam  positions  with  three  digit  vectors,  use  the  Section[54|basis- vector  tiling 
twice  to  create,  in  the  example  here  with  Ri„  =  [  J  ~2  ]  and  Rout  =  [')  , 

containing 

an  anchor  and  points 

point  from  exactly  from 

small  tiles'^  intermediate  lattice  ^2R0Ut1B+  ^  |Ri„|  ^  steering  lattice  Z2Rh[1R~t1B+  and 
large  tiles  ^  array-factor  periodicity  lattice  Z2B+  ^  Roul  |  ^  intermediate  lattice  S2R“11B+. 
Each  beam  position  can  then  be  uniquely  expressed  as 
a  large- tile  anchor  point  *B+ 

+  an  offset  [j]Rout1B+  to  one  of  that  large  tile’s  small-tile  anchor  points 

+  an  offset  [fe]Ri^1RtM1B+  to  one  of  that  small  tile’s  beam  positions. 

The  (example)  has  i  =  [—1  —1],  =  [12  8],  and  [fc]  =  [1  1],  Beam  outputs  are  taken  from  FFT 

output  bins  in  theffiiglTligKtecr^gionj  located  by  setting  i  =  0  to  address  bins  with  [j]  and  [fc]  only. 
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where,  if  we  use  the  optional  but  convenient  basis-vector-tile  approach  of  Scction[53]abovc  (which  approach 
is  itself  based  on  Section  [533]), 


k'=  (60) 

\M  =  k-  k' Rin. 

The  decomposition  can  be  determined  graphically,  using  a  diagram  like  that  of  Fig.  [IT]  or  an  approach 
analogous  to  robust  anchor-point  index  computation  (|4T|)  can  be  used  to  compute  small-tile  anchor-point 
index  ([60])  equally  robustly.  Either  way,  the  values  that  can  be  taken  by  E  are  the  |  Rin  |  specific  integer  row 
two-vectors  for  which  ER^ 1  [  =  0. 

Second,  the  basis  vectors  of  array-factor  periodicity  lattice  Z2B  :  can  be  used  in  a  similar  way  to  con¬ 
struct  large  beamspace  tiles,  as  pictured  in  the  example  of  Fig.  m  to  tile  the  plane  such  that  each  tile  contains 
one  point,  its  anchor  point,  from  that  array-factor  periodicity  lattice  and  such  that  each  point  of  intermediate 
lattice  Z2R(;u!  B  "  is  in  exactly  one  tile.  This  leads  to  a  unique  decomposition  of  an  arbitrary  point  from  that 
intermediate  lattice  as  the  position  of  an  anchor  point  plus  an  intra-tile  offset  in  either  of  the  two  equivalent 
forms 


*'R^B+  =  iB+  +  @RjB+ 
k'  =  iRout  +  [2 

where,  again  using  the  convenient  basis-vector-tile  approach  of  Section [53] 

i  =  Lfc'RouiJ , 

\j\  =  k'  -  4 Rout- 


(61) 

(62) 


(63) 


The  decomposition  can  be  determined  graphically,  using  a  diagram  like  that  of  Fig.  [IT]  or  an  approach 
analogous  to  robust  anchor-point  index  computation  (|4T|)  can  be  used  to  compute  large-tile  anchor-point 
index  ([63])  equally  robustly.  Either  way,  the  values  that  can  be  taken  by  [j]  are  the  |Rout|  specific  integer  row 
two-vectors  for  which  [[jjRoutJ  =  0- 


Small-tile  decompositions  ([58])  and  (59)  can  be  combined  with  large-tile  decompositions  (|6T|)  and  (62) 
to  yield  the  decompositions  of  the  beam  position  and  the  beam-position  index  vector  given  by 


fcRr1R0U]B+  =  ,B  +  0R-|B++  ER^R(7ulB+  (64) 

k  =  4  Rout  Rin  +  [j]Rin  +  El  ■  (65) 


In  Fig.[TT]the  three  thin  line  segments  leading  from  the  origin  to  the  one  (■example  )  beam  position  correspond 
to  the  three  terms  on  the  right  in  beam-position  decomposition  ([64[i-  (The  fourth  segment  is  discussed  below.) 


A  beam  position  is  indexed  or  addressed  by  row  vector  k  or,  using  address  decomposition  ([65]),  the  three 
row  vectors  i,  \j\,  and  E,  so  in  effect  ([65])  relates  two  beam-position  addressing  schemes.  It  also  generalizes 
on  mixed-radix  representation  of  an  integer.  Here  i,  [j],  and  E  are  digit  vectors  from  Z2and  can  respectively 
take  an  infinite  number  of  values,  exactly  |Rout|  values,  and  exactly  |R;n|  values.  Resampling  matrices  Rm 
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and  Rout  are  radix  matrices.  The  most  significant  vector  (MSV)  is  i,  and  the  least  significant  vector  (LSV) 
is  \M-  When  the  vectors  and  matrices  involved  are  replaced  with  integers,  the  procedure  here  for  obtaining 
the  digit  vectors  in  fact  degenerates  to  the  usual  procedure  for  obtaining  a  mixed-radix  representation  of  an 
integer. 

The  key  ideas  of  this  beamspace  tiling  are  summarized  for  review  convenience  in  the  Fig.[TT|caption. 
5.6.5  Two-Tiling  Element  and  Beam  Addresses  Factor  a  Generalized  DFT 

The  General  Form 

Decompositions  ([57])  and  (|65|)  respectively  express  the  element-position  and  beam-position  index  vectors 
as  the  three  terms  above  and  left  of  the  double  lines  in  Table[2]  making  the  product  fcR^c'ern  =  fcR^ 1  Rou[  n 
just  the  sum  of  the  nine  terms  on  the  lower  right  of  the  double  lines.  The  highlighted  terms  are  integers  and 
so  do  not  contribute  to  the  complex  exponential 


e— j27rfcR ste^ra  _  e-j27r|fc]Rin1|m|  g-j27r|fe]Rin1Rout1@  e-j27r[T]R0Ut1[n| 

T>  1 

.steer 

Substituting  this  complex  exponential,  element-position  decomposition  ([57]),  and  beam-position  decom¬ 
position  ([65])  into  beam  computation  ([28])  yields,  after  some  reordering  and  three  quite  optional  uses  of 
Rsteer  =  RoutRin  that  simply  affect  how  the  result  is  written, 


s° 

beam  i Rsteer +[j]R-m  +  IS 


(66) 


|  Rill  I  generalized  |Roul|  generalized 

output  DFTs  input  DFTs 

(one  for  each  [ ~k\)  twiddle  factor  (one  for  each  [77] ) 

Y  A  V  /p-j27r®R7eerEl  ~  j  27rlR“  ^ 

A0  E1  =  A  e  A  •CI53,I3 e 

[n]  such  that  y  [m]  such  that 

LR„ZsJ=o  IR^UnlM 


e-j27r[j]Rollt1[n|^ 


(67) 


fl'lrol.lnl 


A 


£[*-»»°»] 


n  =  Rstcer^+RoullZZl  +  lll 


(68) 


Table  2  —  Terms  of  fcRste'ern  When  Two  Tilings  are  Used  in  Each  Space 


(one  of  DR;,,1  Rout  (one  of  -►) 

-f^out-E^in^ 

Rout  l^'l 

M 

^-R'out-f^in 

^rEoutRjn£ 

*  Rout  (Ml 

iM 

[7]Rin 

\j\nm£ 

\j\m 

SRoul® 

M 

M£ 

[Er^'m 

Er^1rou!ei 
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Here  ([68])  in  effect  addresses  elements  using  two  column-vector  digits  to  generalize  relationship 
governing  the  folding  (data  turning)  of  weighted  element  outputs  into  FFT  input  bins.  Very  general  FFT 
factorization  ([67])  then  transforms  the  doubly  indexed  input-bin  signals  ^  of  (68 1  to  create  doubly  in¬ 
dexed  output-bin  signals  X ^  g.  The  latter  formally  become  beam  outputs  in  ([66]),  which  is  a  generalization 
of  (42 1  to  address  beams  with  two  row-vector  digits. 


It  is  natural  to  suppose  that  the  choice  of  tiling  strategies  somehow  is  key  here,  but  it  isn’t  actually  so. 
Factorization  Rsteer  =  RoutRin  implies  dual  intermediate  lattices  in  position  space  and  beamspace,  and  the 
selection  of  that  dual  pair  of  lattices  indirectly  designs  FFT  factorization  (f67]>.  Choosing  tilings  helps  us  get 
from  the  intermediate  lattice  to  this  factorization,  but  the  factorization  that  results  doesn’t  actually  depend 
on  the  specific  tilings  chosen.  Using  basis-vector  tilings  led  to  expressions  [RoutElJ  =  0  and  |R~ 1  [m]J  =  0 
in  FFT  ([67]),  but  that  tiling  choice  was  purely  for  convenience.  It  gave  us  a  simple  way  to  compute  and  label 
indices.  Prototile  references  more  general  in  form  could  have  been  used,  however,  without  changing  the 
resulting  FFT  structure  in  any  way,  because  the  specific  tilings  chosen  actually  affect  bin  labeling  only.  It  is 
the  intermediate  lattice  used  that  determines  the  FFT  factorization  itself. 


The  Block  Diagram 


FFT  expression  ([67])  is  presented  as  a 


3D  block  diagram  in  Fig.[l2]for  the  example  case  in  which 


j  input  msv  m\  e  {o,  m\  m2}  c  % 2 

\  output  LSV  \k\  G  {0,  H1,  [fc]J  }  C  5E2, 


r  input  lsv  m  e  {o,  m\ m\  @3}  c  %2, 

\ output  MSV  \j\  G  {0,  \j\\  [j]2,  \j\3}  C  Z2 


This  diagram  is  very  similar  in  appearance  to  Fig.  [9]  and  the 
a  large  DFT  in  terms  of  smaller  ones.  In  other  ways  the  two 
considered  in  its  full  generality. 


two  systems  have  in  common  that  they  express 
systems  are  quite  different,  at  least  if  Fig.  12  is 


Fig.  12  —  General  FFT  digital  steering  structure,  here  with 
|Rin|  =  3  and  |  Rout  |  =  4,  transforms  weighted  and  folded 
element  outputs  into  beam  outputs  A'jjj  gj  in  steps: 

1.  a  ruralized  DFT  for  each  [n]  transforms 
position  index  [m]  to  beamspace  index  0  using  R„, 

2.  a  twiddle-factor  multiply  per  0,  El  combination 
using  Rj„  and  Roui,  a  no-op  when  [0  =  0  or  E  =  0,  and 

3.  a  generalized  output  DFT  for  each  [0  transforms 
position  index  E  to  beamspace  index  [J]  using  Rout- 


In  Fig.  [9]  all  component  transforms  were  ID  FFTs  with  integer-indexed  bins,  but  here  in  Fig.  [12] each 
component  transform  is  a  generalized  DFT  with  inputs  and  outputs  indexed  by  vectors  from  ff1  and 


respectively.  In  the  example  application  addressed  by  the  maps  of  position  space  and  beamspace  in  Figs.  10 
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and|TTj  index  vectors  \m ]  and  \]\  select  small  tiles  related  to  an  intermediate  lattice,  and  index  vectors  \n\  and 
| fcj  select  bins  from  inside  those  tiles.  This  is  very  different  in  principle  from  the  row/column  indexing  of 

Fig-El 


Here  the  component  transforms’  vertical  and  horizontal  orientations  in  the  diagram  and  the  arrangement 
of  input  and  output  bins  on  grids  need  bear  no  particular  relation  to  the  geometrical  positioning  of  input 


and  output  bins  in  position  space  and  beamspace.  In  Fig.  12  the  row-column  arrangement  of  inputs  and 
outputs  is  largely  just  a  drawing  convenience  to  help  make  the  factorization  comprehensible,  though  there 
is  a  pattern  to  it:  for  both  input  and  output  transforms  changing  the  least  significant  digit  vector  steps  from 
one  component  transform  to  another,  while  changing  the  more  significant  digit  vector  steps  from  bin  to  bin 
within  a  single  component  transform. 


12 


multiplications  by  twiddle  factors  — this  is  actually  the  traditional 


Finally,  in  FFT  (|67j)  and  in  Fig. 
name  in  the  FFT  literature — come  between  the  two  layers  of  smaller  DFTs.  This  twiddle  factor  is  not  shown 


when  =  0  or  \n\  =  0  because  in  those  cases  its  value  e 


-j27r[fc]Rst 


=  1. 


5.7  Three  Key  Special  Cases 


Technically  the  DFT  factorization  into  smaller  DFTs  of  Fig.  12  can  only  be  called  an  FFT  if  each  of 
its  component  DFTs  is  either  itself  an  FFT,  of  one  kind  or  another,  or  a  very  small  DFT,  small  enough  that 
using  an  FFT  approach  to  calculating  it  would  be  unwarranted.  Three  specializations  of  particular  interest 
correspond  to  three  specific  factorizations. 


5. 7. 1  The  Simple  Case 


The  simple  case  of  Section  [575]  above  results  from  substitutions 


H0i  i:  — 


1  0 

0  iVe? 

rri\  ~ 
0 

0 

lm,2  J 


@  = 

Rin  = 

E  = 


[0  r2] , 

Xx  0 

0  1 

[ri  0] 


where  each  of  r\,  r2,  mi,  and  m2  range  over  the  integers  0, . . . ,  Nex  —  1.  These  turn  the  general  FFT 
realization  presented  here  in  ([67])  and  Fig.  [12]  into  Section  [53] s  special-case  FFT  realization  of  ([45])  and 
Fig.  [9]  where  input  and  output  FFTs  are  just  ordinary  ID  FFTs  of  Nex  points  each. 


5. 7.2  Coarse  and  Fine  Beam-Spacing  Modes 


Coarse  and  fine  beam-spacing  modes  can  be  realized  using  an  extra  small-determinant  factor  as  the 
rightmost  factor  in  Rsteer  as  in  the  example  of  (|48[).  These  substitutions  make  Rin  such  an  “extra”  factor 
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with  Rout  the  “main  portion”  of  Rst 


Rout  — 


iVex 

0 

0  ' 

iVex. 

9 

0 

=  [r-1  r2\ 

0,  m\--- 

Rin 

small, 

mi ' 
m2 

5 

E 

g  {o,E\. 

This  gives  the  general  FFT  structure  of  Fig.  12  a  large  number 
small  number 


I  Rout  |  =  iVgX  of  small  inpu  FFTs  and  a 


Rin  |  of  large  output  FFTs.  Each  of  the  latter  large  FFTs  are  implemented  as  in  Fig.  [9]  or, 
equivalently,  using  the  simple  case  described  above.  Each  of  n,  r2,  mi,  and  m2  here  again  range  over  the 
integers  0, . . . ,  Nex—  1,  and  the  mismatch  here  between  these  variable  names  and  the  names  of  the  vectors 
they  make  up  is  simply  for  compatibility  with  Fig.  [9] 

Here  digit  vectors  \fn\  and  E  each  take  Rm  specific  values,  one  of  which  is  a  zero  vector.  The  other 


specific  values  can  be  determined  graphically,  as  illustrated  in  Figs.  10  and  11  for  |Rjn|  =  3,  or  a  some 


what  larger  set  determined  graphically  can  be  thinned  down  computationally  to  keep  only  those  ]m]  and  E 
for  which  [R~uj  [mlj  =  0  and  ERjn  1  =  0  respectively,  with  these  floor  tests  computed  using  the  robust 


approach  of  ([35])  and  (|4T])  that  was  derived  above  in  Section  5.3.4 


The  example  design  presented  Figs.  [K4  and  11  for  example,  would  be  realized  using  Fig.  12  with 


Hr  =  256  three-input,  three-output  generalized  DFT  blocks  as  inpul  FFTs  and  three  16  xl6  point  2D  FFTs 
as  output  FFTs.  Each  of  those  three  2D  output  FFTs  would  be  realized  using  32  ID  FFTs  of  16  points  each 
using  the  architecture  of  Fig.  [9]  The  Appendix  shows  that  the  three-point  generalized  DFTs  are  actually 
ordinary  ID  DFT  (or  FFT)  blocks. 


5.7.3  Efficient  Computation  of  Beam  Clusters 


Efficient  computation  of  beam  clusters  is  often  possible  using  either  the  first,  simple  case  above  or  using 
an  Rsteer  =  RoutRin  factorization  in  which  Rout  is  the  small  “extra”  factor,  basically  reversing  the  roles  of 
Rin  and  Rout  relative  to  the  case  just  discussed.  Now 


Rout 

small, 

\m\  = 

'mj ' 

m2 

m  g  jo,®1, . . . , @|Roi,t|  1}, 


H  efoj1 . i11™*'-1}, 

'iVex  0 


Rin  = 


vex 

0  AH 


E  = 


[ri  r2] . 


The  difference  between  this  case  and  the  last  is  one  of  input-bin  and  output-bin  geometry  in  position  space 
and  beamspace.  The  last  case  corresponded  to  Figs.  [10]  and  [TT]  in  which  input  bins  were  structured  into  a 
few  giant  blocks  and  output  bins  were  structured  into  many  small  ones.  Here  that  is  reversed  so  that  input 
bins  are  structured  into  many  small  blocks  and  output  bins  into  a  few  large  blocks.  This  is  convenient  for 
computing  a  cluster  of  beams,  because  beams  near  each  other  in  beamspace,  such  as  those  in  a  cluster, 
are  often  computed  near  each  other  in  the  same  large  block  of  output  bins.  This  makes  it  easy  to  improve 
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computational  efficiency  by  deleting  entire  unused  rows  or  columns  of  output  bins  in  the  large-block  FFTs. 
Examples  of  such  systems  are  presented  in  Ref.  [3] 


6.  SUMMARY  AND  OBSERVATIONS 

This  document  used  a  simple  theoretical  formulation  of  a  planar  array  of  identical  elements  positioned 
on  points  of  a  lattice  and  surrounded  by  guard  elements  to  develop  many-beam  phase-shift  steering  using  a 
generalized  Cooley-Tukey  FFT  structure  to  realize  a  generalized  2D  DFT,  one  in  which  the  usual  number  of 
DFT  points  has  been  replaced  with  a  factorable  2x2  resampling  matrix  Rsteer- 

In  discussions  of  FFT  beamsteering,  four  specific  questions  typically  arise,  two  that  are  about  the  FFT 
itself  and  two  that  concern  its  application  to  radar  systems. 


1.  Are  these  “decimation  in  time”  or  “decimation  infrequency”  FFT  structures? 

Such  names  only  apply  here  beginning  in  Section [576|  when  an  intermediate  lattice  is  first  used,  and  even 
then  only  when  one  of  |R;n|  and  |Rout|  is  small  enough  that  the  corresponding  resampling  matrix  is  used 
without  further  factorization,  when  the  associated  DFTs  are  realized  directly  from  general  expression 
((43]).  In  that  case  the  classification  is  according  to  which  of  |Rin|  or  |Rout|  is  large  and  which  is  small.  In 
Table  [3] position  corresponds  to  the  usual  “time’,’  and  beamspace  corresponds  to  the  usual  “frequency.” 


Table  3  —  Decimation  in. . .  What? 


structure 

1  Ri„  | 

|  Rout  I 

examples 

decimation  in  position 

large 

small 

Section 

5.7.3 

Efficient  Computation  of  Beam  Clusters 

decimation  in  beamspace 

small 

large 

Section 

exampk 

5.7.2 

facte 

Coarse  and  Fine  Beam-Spacing  Modes, 
trization  (|48]),  example  tilings  of  Figs.  10  and  11 

2.  What  is  the  purpose  of  the  tiles? 

The  short  version  is  that  large  input  tiles  are  about  input  folding  (data  turning),  large  output  tiles  are 
about  beam-position  periodicity,  and  small  input  and  output  tiles  are  what  enable  the  DFT  to  be  turned 
into  an  FFT  through  factorization.  The  general  tiling  rules  retain  the  full  design  flexibility  inherent  in  the 
DFT  and  FFT  mathematics  while  organizing  matters  graphically  to  make  comprehension  less  arduous. 
The  longer  version  goes  something  like  this. 

Large  input  tiles  structure  the  discrete-space  Fourier  transform  into  a  finite  DFT  sum  and  an  infinite 
input-folding  sum,  in  the  process  giving  the  DFT  a  finite  number  of  input  bins,  by  making  product 
km  an  integer  in  the  derivation  of  the  multibeam  computation  of  ([36])  and  (f37j).  Further,  they  rep¬ 
resent  the  most  general  way  of  making  that  happen.  The  tiles  then  imply  an  MSV-LSV  addressing 
scheme  that  maps  any  element  location  to  a  combination  of  an  DFT  input  bin  and  an  input-folding 
(data  turning)  term. 


54 


Jeffrey  O.  Coleman 


Large  output  tiles  eliminate  redundant  DFT  output  bins  by  making  product  j  [n]  an  integer  in  the  deriva¬ 
tion  of  the  finite  set  of  beams  computed  in  ([42])  and  (|43j).  They  represent  the  most  general  way  of 
making  that  happen.  The  tiles  then  imply  an  MSV-LSV  addressing  scheme  that  maps  any  beam 
position  to  a  combination  of  an  DFT  output  bin  and  a  steering- vector  period. 

Small  tiles  add  “intermediate-significance  vectors”  to  create  MSV-ISV-LSV  element-location  address¬ 
ing  and  MSV-ISV-LSV  beam-position  addressing  in  such  a  way  that  the  DFT  factors  into  LSV- 
labeled  nput  FFT  blocks  with  ISV  input-bin  labeling  followed  by  (twiddle-factor  multiplications 
and)  LSV-labeled  output  FFT  blocks  with  ISV  output-bin  labeling.  Fundamentally,  the  tiling  scheme 
enables  this  factorization  by  making  [7]lml.  the  middle  term  of  the  nine  tabulated  in  Table  [2]  an  in¬ 
teger.  It  is  the  most  general  way  of  making  that  happen.  The  specific  FFT  factorization  realized  is 
determined  by  the  choice  of  the  intermediate  lattice  and  its  dual  used  in  the  construction  of  the  small 
tiles. 

3.  What  are  the  main  limitations  of  the  approach? 

The  key  limitation  of  the  approach  is  that  the  positions  of  the  synthesized  beams  in  beamspace  (sine 
space,  give  or  take  a  scale  factor)  must  be  drawn  from  a  superlattice  of  the  array-factor  periodicity 
lattice,  with  the  two  lattices  related  by  Rsteer-  Further,  for  an  FFT  structure  to  be  viable  Rsteer  must  be 
factorable  into  resampling  matrices  with  determinants  small  in  absolute  value.  The  FFT  incorporates  a 
specific  choice  of  Rsteer  in  a  deep,  structural  way,  so  generally  Rsteer  cannot  be  arbitrarily  changed  in  real 
time  and  therefore  neither  can  the  steering  lattice  of  beam  positions.  At  most  there  might  be  flexibility  to 
switch  between  steering  lattices,  particularly  between  related  ones,  for  example  as  in  the  case  of  “coarse 
and  fine  beam-spacing  modes”  discussed  near  the  end  of  the  last  section.  But  these  restrictions  do  mean 
that  not  even  the  beamspace  spacing  of  the  steering-lattice  beams  can  be  freely  chosen.  Enough  viable 
choices  of  Rsteer  are  available,  however,  that  respecting  these  limits  seldom  cripples  system  design,  as 
long  as  specific  beam  positions  are  less  important  than  having  the  spacing  between  them  be  small  enough 
without  being  too  small. 

4.  When  this  approach  is  used  for  a  shipboard  array,  can  it  be  modified  to  include  electronic  ship-motion 
compensation  ? 

Yes  and  no.  Yes,  it  is  easy  enough  to  replace  element  outputs  s°  with  s[)r.  /Atomp  everywhere,  so  that 
specific  element  output  is  replaced  with  e_j27r^A  °mp”  using  one  extra  complex  multiply  per  element 
to  effectively  translate  the  entire  lattice  of  beam  positions  by  kAmp  =  /AcompB+  to  compensate  for 
motion.  But  no,  this  is  not  complete  compensation.  One  beam,  at  least  if  it  is  circularly  symmetric  in 
beamspace,  can  be  compensated  by  a  simple  translation  in  this  way.  But  translation  is  not  enough  when  a 
family  of  beams  on  lattice  points  is  involved,  because  there  is  no  way  to  compensate  for  ship  rotation,  for 
example  ship  roll  when  a  radar  face’s  boresight  vector  points  forward  or  ship  pitch  when  a  radar  face’s 
boresight  vector  points  abeam.  Even  when  motions  are  small  enough  that  the  trigonometric  distortions 
involved  in  mapping  from  ship-motion  angles  to  changes  in  beamspace  position  k  are  not  significant, 
there  is  a  problem:  rotation  in  beamspace  of  the  entire  set  of  beams  would  be  needed  for  complete 
motion  compensation,  and  this  is  simply  not  possible.  Scatterer  positions  must  simply  be  estimated  in 
ship  coordinates,  converted  to  Earth  coordinates,  and  passed  in  that  form  to  the  tracking  system. 
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A  PRIME  |R|  MAKES  A  GENERALIZED  2D  DFT  INTO  AN  ORDINARY  ID  DFT 
INTRODUCTION 

A  DFT  Configuration  Seen  Throughout  This  Document 

Every  DFT  seen  heretofore  in  this  document,  whether  a  “top  level”  DFT  or  a  component  DFT  inside  an 
FFT  expression,  is  of  the  form 

Xm=  Y.  (Al) 

iG|2 such  that 
Bloc®  e  (input  prototile) 


for  some  basis  matrix  Bioc  and  some  resampling  matrix  R,  where  E>  x^,  and  are  variables  local  to  this 
discussion  and  may  or  may  not  be  the  specific  array-related  values  from  before.  Here  (input  prototile)  is  a 
position-space  tile  containing  the  origin  such  that  tiling  all  of  position  space  with  its  translates  gives  each 
tile  exactly  one  anchor  point  from  sublattice  BiocR^2  and  exactly  |R|  points  from  lattice  B|oc^2.  Fikewise, 
we  have  supposed  E  €  IE2  such  that  EB^e  (output  prototile),  where  (output  prototile)  is  a  beamspace  tile 
containing  the  origin  and  where  tiling  all  of  beamspace  with  its  translates  gives  each  tile  exactly  one  anchor 
point  from  Z2B|oc  and  exactly  |R|  possible  intra-tile  offsets  from  superlattice  Z2R_1Bj{^.  We  often  used 
basis- vector  tiling  in  particular,  which  leads  to  prototile  descriptions 

\n\  such  that  [_  R_1[n]J  =  0  , 

E  such  that  LER'1  J  =  0 

but  those  particular  structures  are  completely  optional. 

The  Result  to  be  Proved 

When  |R|  is  prime,  the  sets  of  allowed  index  vectors  M  and  E  for  the  intra-tile  offsets  can  each  be  put 
in  some  order 

{allowed  [n]}  =  {[n]°  m\  EST •  ■  • ,  (A4) 

{allowed  E}  =  {E°  E{  E{  •  •  • ,  E|R'-1}  (A5) 

such  that,  if  we  use  xn  for  x^n,  use  X k  for  X^k ,  and  let  A  =  |R| ,  the  generalized  2D  DFT  in  (|A1|)  becomes 
the  ordinary  one-dimensional  DFT 

N- 1 

Xk=Y*ne-^. 

n=0 


(A2) 

(A3) 
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This  allows  the  implementer  the  convenience  of  drawing  on  canned  FFT  routines,  since  for  a  prime  number 
of  points  the  FFT  is  just  a  straightforward  DFT  sum. 

The  proof  offered  below  is  constructive,  proving  not  only  that  such  orderings  exist  but  providing  a 
method  for  determining  them.  In  fact  |R|  —  1  choices  of  these  orderings  are  provided. 

A  BACKGROUND  FACT 

The  proof  below  hinges  on  this  result  from  elementary  number  theory. 

Bezout’s  identity.  If  k\  and  k2  are  nonzero  integers  with  gcd(fci,  k2 )  their  greatest  common  divisor, 
there  exist  integers  n\  and  n 2  such  that 


kin\  +  k2n2  =  gcd(Aii,  k2). 


One  special  case  is  of  particular  interest.  If  integer  M>  1  with  M  prime,  then  gcd(fc,  M)  =  1  for  any  k  that 
is  not  a  multiple  of  M  or,  equivalently,  for  which  k  mod  M/0.  Bezout’s  identity  then  tells  us  that  we  can 
find  integers  —i  and  n  such  that  these  three  equivalent  statements  hold. 


kn  —  iM  =  1 , 
kn  =  £M  +  1, 
kn  mod  M  =  1. 


The  latter  reveals  that  Bezout’s  identity  is  telling  us  that  every  k  that  is  nonzero  modulo  M  has  a  multiplica¬ 
tive  inverse  modulo  M.  Of  course  this  also  follows  from  another  algebraic  fact,  that  the  integers  modulo  a 
prime  form  a  field. 

Multiplicative  inverses  or,  more  generally,  Bezout  coefficients  can  always  be  found  with  the  extended 
Euclidean  algorithm  Q,  but  when  the  solutions  sought  are  restricted  to  the  modulo  M  integers  with  M 
small,  as  is  generally  the  case  in  this  application,  a  brute-force  computational  search — just  try  all  possible 
solutions — is  a  simpler  and  quite  reasonable  alternative. 

THE  PROOF 

First:  Express  the  Elements  of  R,  _ 1  as  Rational  Numbers 

Recall  that  the  inverse  of  a  nonsingular  matrix  R  can  be  expressed  in  terms  of  its  adjugate  matrix  and 
its  determinant  as 


(A6) 


Of  course  any  number  is  the  product  of  its  absolute  value  and  sign,  so  det(R)  =  |R|  sgn(det(R))  = 
|R|  /  sgn(det(R)).  Combining  with  (|A6[)  allows  us  to  write 


(A7) 


Rj  =  sgn(det(R))  adj(R). 
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Resampling  matrix  R  is  an  integer  matrix  and  therefore  so  is  R, . 

Second:  Tile  Beamspace,  Pick  a  Seed  Point  in  Beamspace ,  and  Analyze  That  Point ’s  Coordinates 

Lattice  Z2Bjj„  has  a  superlattice  Z2R_ 1  B|nc.  Partition  the  plane  into  identical  tiles,  each  containing 
one  anchor  point  of  lattice  Z2Bj^c.  Then  each  tile  contains  exactly  |R|  points  of  superlattice  Z2R'  1  B]oc. 
Choose  any  “seed”  point  fcseedR“1B]J)c  in  Z2R~ 1  B](l)c  that  is  not  in  Z2B^c.  Of  course  there  arc  |R|  —  1 
possible  choices  of  seed  point  in  each  tile,  and  this  is  ultimately  why  we  end  up  with  |R|  —  1  choices  of 
orderings  (|A4[)  and  (|A5[).  Using  this  new  seed,  let 

fcmod  =  /bseedR.  mod  |R|5  (A8) 

with  the  modulo  operation  taken  on  each  element  of  the  row  vector  fcseedRj  separately.  That  definition  of 
kmod  is  equivalent  to  the  first  of  these,  and  the  second  then  follows  from  (|A7[): 

fcseedR.  =  j|R|  +  kmod  for  some  j  €  Z2,  (A9) 

fcseed R-i  =j  +  _Lfc™d  (for  that  same  j). 

|R| 

In  the  latter  we  express  the  elements  of  row  vector  fcseedR_1  as  the  sum  of  their  integer  and  fractional  parts. 
All-integer  equivalent  (|A9[)  is  easier  to  work  with  below,  however.  While  we  could  easily  compute  row 
vector  j,  actually  doing  so  is  unnecessary.  It’s  enough  to  know  that  such  a  vector  exists. 

One  useful  fact  about  A:mod  is  easily  established.  Since  kseedR~  1  Bj£c  is  not  in  Z2B^)c,  equivalent 
statements 


fcseedR-1B+c/jB+, 

feseedR— 1  /j, 

kseedRi  |R| 


all  hold.  Comparing  the  last  of  these  to  (A9i  reveals  that  fcmod  /  [0  0]. 

Third:  Construct  a  Seed  Point  in  Position  Space  with  a  Required  Property 

We  can  now  construct  a  nonzero  seed  point  Biocnseed  in  position  space  in  such  a  way  that  these  equiva¬ 
lent  statements  are  true. 


fcseed R,  nseed  =  ( j  nseed  +  £)  \  R|  +  1  for  some  integer  t, 
fcseed =  j  nseed  +  l  +  ,  1  (for  that  same  (j. 

R 


(A  10) 
(All) 


The  second  statement  decomposes  the  scalar  fcscedR,  1nseed  into  integer  and  fractional  parts.  What  is  im¬ 
portant  to  us  is  the  particular  value  1/|R|  of  that  fractional  part.  It  is  easiest  in  establishing  this,  however, 


to  aim  for  (A10),  the  all-integer  statement  of  the  same  result.  Consider  two  distinct  cases,  according  to  the 


form  of  kmod  from  dX8b. 
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Case  I:  fcmod  =  [k\  k2]  with  exactly  one  of  k  \  and  k2  equal  to  zero. 

Suppose  one  of  k\  and  k2  has  value  0  and  the  other  has  value  k  with  0  <  k  <  |R|.  Let  n  with 
0  <  n  <  |R|  be  the  multiplicative  inverse  of  k  modulo  |R|  so  that  kn  =  1  +  f'|R  for  some  integer  £. 
Let  nseed  =  [  o  ]  or  nseed  =  [  °  ]  according  as  fcmod  =  [k  0]  or  femod  =  [0  /,:] .  Using  (|A9|)  establishes  (] A 1 0[) : 


fcseedR;nseed  =  jnseed|R|  +  fcmodnseed 
=  jnseed|R|  +  kn 
=  jnseed|R|  +  (1  +  £|R|). 


Case  II:  fc 1,10(1  =  [k\  k2]  with  both  k\  and  k2  nonzero. 

Special  case:  gcd(fci,  k2)  =  1.  This  is  a  special  case  of  the  general  case  considered  next,  and  its  proof 
is  included  only  as  an  easier-to-follow  preview  of  the  proof  of  the  general  case  to  follow. 

By  Bezout’s  identity  there  are  integers  n±  and  ri2,  readily  found  in  fact,  such  that  k )  n  \  +  kyria  =  1. 
Let  nseed  =[”’].  Then,  by  @5), 


fcseedRinseed  =  j  nseed|R|  +  hm  +  k2n2 


establishing  (|A1Q[)  with  £  =  0. 


General  case:  Let 


K  =  fei/gcd(fei,  k2), 
k2  =  k2/gcd(ki,k2) 


so  that  integers  k[  and  k2  have  gcd(fc/1,  k'2)  =  1  and 

k\  =  k[  gcd(ki,k2), 
k2  =  k'2gcd(ki,  k2). 

Let  ngcd  be  the  multiplicative  inverse  of  gcd(&i,  k2)  modulo  |R|.  Then  there  is  some  integer  l  such 
that  gcd(/ci,  k2)ngcd  =  1  +  ^|R|,  which  with  the  above  gives 

fcingcd  =  k[  gcd(/ci,  k2)ngcd  =  k[(  1  +  £|R|), 
k2ngcd  =  k2gcd(ki,k2)ngcd  =  A4(1  +  ^|R|). 


By  Bezout’s  identity  there  are  integers  n\  and  n2,  readily  found  in  fact,  such  that  k\  n  \  +  k2n2  =  1. 
Let  nseed  =  [  ] .  By  (|A9[)  and  then  the  above, 

fcseedRinseed  =  jnseed|R|  +  hngcAn  1  +  k2ngcdn2 

=  jnseed|R|  +  (k[m  +  k2n2)(l  +  €|R|) 


establishing  (|A1Q  >. 
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Fourth:  Construct  Invertible  Functions  from  the  Modulo  |R|  Integers  to  the  Vectors  in  the  Prototype  Tiles 
For  each  n  E  {0, . . . ,  |R|  —  1}  find  the  anchor  point  B]0CRraand  the  intra-tile  offset  B|oc[n]  such  that 

nBiocnseed  =  Bi0CRm+  Bioc[rg,  (A  12) 

mrseed  =  Rm+  \n\  (A  13) 


where  the  second  is  of  course  obtained  by  the  first  by  premultiplying  by  B  and  using  identity  Bj^cBi0C  = 
I.  Suppose  the  value  of  function  fin)  is  given  by  \n\  of  the  above  decomposition. 

Is  function  fin)  invertible?  Given  some  allowed  [n]  and  relationship  (bin)  =  [to],  is  n  established  unam¬ 
biguously?  It  is  if  we  can  show  that  fin1)  =  f(n2)  implies  n1  =  n2.  Suppose  fin1)  =  \n\  =  fin2)  so 
that 


n 1nseed  =  Km1+m: 
n2  nseed  =  Rm2+  M- 


Subtract  the  two  equations  to  obtain,  in  three  equivalent  forms, 

(n1  —  n2)  nseed  =  R  (m1-™2), 
(n1-n2)R“1nseed=  (m1-™2), 
(n1— n2)Rj  nseed  =  iRKm1— m2). 


By  the  last  of  these  |R|  is  a  factor  of  both  elements  of  vector  (n1  —  n2)R(nseed,  which  means,  since  |R|  is 
prime  and  cannot  be  split,  it  must  either  be  a  factor  of  integer  n1  —  n2  or  of  both  elements  of  vector  R,nseed. 
But  the  latter  would  make  integers  of  both  elements  of  JjtR(  TOseed  =  R 


- 1 ^  seed 

'  V  1 


and  this  would  in  turn 


require  fcseedR.  1  nseed  to  be  an  integer  and  so  contradict  (|A1 1  [).  Therefore  |R|  divides  integer  n1  —  n2.  But 
n1and  n2are  each  restricted  to  0, ,  |R|  —  1,  so  n1  —  n2  E  {— (|R|  — 1), . . . ,  |R|  —1}.  The  only  multiple 
of  |R|  in  the  latter  set  is  zero,  so  n1  —  n2  =  0  or  n1  =  n2. 


The  coiTespondence  between  n  and  position-space  vector  [to]  just  established  has  a  counterpart  in  a 
relationship  between  k  and  beamspace  vector  \M-  Suppose  that  for  each  k  E  {0, . . . ,  |R|  —  1}  anchor  point 
£B^c  and  intra-tile  offset  [fc]R_ 1  B^c  are  the  unique  choices  such  that  these  equivalent  relationships  hold: 


kfcseedR-1B+c=  *B+  +  m-^,  (A  14) 

kkseed  =m  +1.  (A  15) 


Suppose  the  value  of  function  rf{k)  is  given  by  \M  of  the  above  decomposition.  By  an  argument  essentially 
identical  to  that  used  above  for  fin),  the  function  ’fik)  is  invertible.  Given  some  allowed  \M  and  relationship 
fik)  =  integer  k  is  established  unambiguously. 
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Fifth:  Assemble  the  Final  Result 

Given  an  allowed  [n]  there  is  exactly  one  integer  n  with  0  <  n  <  |R|  such  that  f(n)  =  [ra],  so  the  sum  in 
((AT)  can  just  as  well  be  taken  over  those  n.  The  same  terms  are  summed.  Indeed,  we  may  as  well  let  xn 
refer  to  x ^  using  this  same  correspondence.  Similarly,  given  an  allowed  \M  there  is  exactly  one  integer  k 
with  0  <k<  |R|  such  that  f(k)  =  M,  so  we  can  let  Xf.  refer  to  Xjg  using  that  relationship.  Making  these 
substitutions  in  (|A1[)  along  with  substitutions 


M  =  nseedn  -  Rm" 

M  =  kkseed  -  £kR 


from  (|A13[)  and  (|A  15),  where  mand  £  have  been  replaced  with  mn  and  £k  respectively  as  a  reminder  of 


their  dependency  on  n  and  k. 


|R|-i 


_1 

y -j  27rfcfcseed  R- 1  nscedn 


Xk  —  ^  ^  X' 


(A  16) 


n= 0 


But  by  ((ATT), 


]fesecdR-1nseed  =  M  + 


for  some  integer  M,  so  the  remaining  complex  exponential  in  ( A16 1 


so  (A16)  itself  becomes 


|R|— 1 


n= 0 


a  completely  conventional  DFT  of  |R|  points  in  one  dimension. 


Q.E.D. 


A  SIMPLE  EXAMPLE 


The  key  step  was  the  third,  finding  a  suitable  nsecd  to  pair  with  the  chosen  fcseed  so  that  fcseedR,nseed  mod 
|R|  =  1.  Having  established  that  such  an  nseed  exists,  we  can  actually  find  it  in  practice  with  a  simple  brute- 
force  search.  Consider  a  Matlab  example  with 


»  R  =[2,-1; 
»  1,  3] 


>> 


Now  let  us  walk  through  the  four  steps  of  the  proof  but  computing  only  what  is  essential. 
1.  Express  the  elements  of  Rr  1  as  rational  numbers. 
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>>  adR=abs (det (R) ) 
adR  = 

7 

>>  Ri=round ( inv (R) *adR) 

Ri  = 

3  1 

-1  2 

>>  R*Ri  %  a  check 
ans  = 

7  0 

0  7 

2.  Tile  beamspace,  pick  a  seed  point  in  beamspace,  and  analyze  that  point’s  coordinates. 

Begin  either  by  choosing  a  beamspace  seed  point  feseedR”1B1J)C  that  is  not  in  Z2B[jj  or,  equivalently, 
choosing  integer  vector  fcseed  with  fcseedR,  mod  |R|  ^  [0  0] .  We  can  just  guess  and  test  the  result. 

>>  kseed= [-1, -1 ] ; 

>>  mod (kseed*Ri, adR) 
ans  = 

5  4 

3.  Construct  a  seed  point  in  position  space  with  a  required  property. 

Choose  an  nseed  to  make  fcseedRjnseed  mod  |R|  =  1.  Because  this  test  is  done  modulo  |R|,  integer 
multiples  of  |R|  can  be  added  to  the  elements  of  nseed  at  will  without  changing  the  result  of  the  test. 
This  means  it’s  enough  to  just  test  nseed  =  ]  for  each  m,  ri2  G  {0, . . . ,  |R|  —  1}.  (Or  use  any  other  set 

of  integers  no  two  of  which  differ  by  a  multiple  of  |R|.) 

>>  range=  ( 0 :  (adR-1 ) ) ; 

>>  nl=range; 

>>  n2=range; 

>>  (i, j ] =f ind (mod ( repmat (kseed*Ri ( : , 1) *nl . ' , 1, adR) + . . . 

>>  repmat (kseed*Ri ( : , 2) *n2, adR, 1 ) . . . 

>>  , adR) ==1) ; 

>>  nseedlnit= [nl (i ) ; n2 ( j ) ] 
nseedlnit  = 

3  5  0  2  4  6  1 

0  1  2  3  4  5  6 

Of  course  nseed  here  really  represents  position-space  point  Biocnseed,  and  only  this  point’s  intra-tile 
offset  actually  matters.  There  are  |R|  points  in  a  tile,  but  this  calculation  tests  |R| 2  choices  of  ns  ,  so  it 
should  be  no  great  surprise  that  this  test  returns  |R|  suitable  choices  of  nseed. 

If  each  column  output  above  is  taken  to  be  a  choice  of  nseedlml,  it  is  easily  verified  that  each  of  the 
associated  points  B|ocnseedlmt  has  the  same  intra-tile  offset. 

>>  (R*mod (Ri*nseedlnit , adR) ) /adR 
ans  = 

0  0  0  0  0  0  0 

2  2  2  2  2  2  2 

If  we  take  nseedto  be  one  of  the  columns  of  nseedlnit  plus  a  vector  |R|  [  ^  ]>  then  choosing  the  column 
and  mi  and  rn-2  chooses  the  tile.  The  results  do  not  depend  on  which  tile  is  chosen,  so  this  is  only  a 
matter  of  plotting  convenience  (as  in  Fig.|Al|below).  The  chosen  nseed  should  be  tested  to  be  sure. 

>>  nseed=nseedlnit ( : , 2 ) +adR* [ -1 ; 0 ] 
nseed  = 

-2 

1 
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Fig.  A1  —  Graphical  approach  to  ordering  intra-tile  offsets,  which  correspond  to  FFT  input  (left) 
and  output  (right)  bins,  so  that  a  2D  DFT  with  |R|  prime  becomes  an  ordinary  ID  DFT.  Flighlighting 
flags  whether  numbers  indicate  ID  bin  indexes  or  ID  bin  indexes  less  |R|  .  Flere  |Rj  =  7, 
corresponding  to  seven  points  per  tile. 


>>  mod (kseed*Ri*nseed,  adR) 
ans  = 

1 

4.  Construct  invertible  functions  from  the  modulo-  |R|  integers  to  the  vectors  in  the  prototype  tiles. 


According  to  ( A12 1  and  ( A14 1  we  can  obtain  the  ordered  [n]  and  \k\  in  our  |R|  =  7  example  from  the  intra¬ 
tile  offsets  of  B]ocnseedn  and  A:A:seedR_1B]J)C  with  n  and  k  each  taking  values  in  turn  of  0, 1,  2,  3, 4, 5, 6. 
If  we  assume  basis-vector  tiling  has  been  used  in  both  position  space  and  beamspace,  we  can  easily 
compute  the  index  vectors  for  the  intra-tile  offsets. 


>>  nBin= (R*mod (Ri*nseed*range,  adR)  ) / adR 
nBin  = 

0  0  1  1  0  0  1 
0  2  13  13  2 

>>  kBin= (mod ( range . ' *kseed*Ri,  adR)  )  *R/ adR 
kBin  = 

0  0 

2  1 

1  0 

1  2 

2  0 

2  2 

1  1 

Did  it  work?  We  can  just  test  all  the  fch'nR)  rib"1  products  and  see  if  they  are  the  same,  as  hoped,  as  the 
kn  products. 


>>  kbinRinbin=mod (kBin*Ri*nBin,  adR) 
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kbinRinbin  = 


0  0  0 

0  12 

0  2  4 

0  3  6 

0  4  1 

0  5  3 

0  6  5 


0  0  0  0 

3  4  5  6 

6  13  5 

2  5  14 

5  2  6  3 

16  4  2 

4  3  2  1 


>>  all (all (kbinRinbin==mod (range . ' *range, adR) ) ) 


ans  = 

1 


One  other  tiling  scheme  that  is  nearly  as  simple  uses  basis-vector  tiles  offset  by  a  vector.  For  example, 
on  the  right  in  Fig. 


A1 


the  beamspace  tiles  are  basis- vector  tiles  translated  by  — 
vectors  for  the  intra-tile  offsets  is  simple,  because  no  points  lie  on  tile  boundaries: 


loc' 


Computing  index 


>>  offs et=ones( length  (range ) , 2 ) * . 5; 

>>  kBin= round ( (mod (range . '  *kseed/R+of f set , 1) -offset) *R) 
kBin  = 

0  0 

-1  -1 

1  0 

0  -1 

0  1 

-1  0 

1  1 


A  geometric  view  of  the  correspondence  the  above  calculations  establish  between  scalar  indices  n  and  k 
and  the  position-space  and  beamspace  bins  they  are  discovered  to  represent  is  especially  illuminating.  One 
could  plot  left  sides  nBiocnseed  and  kfcseedR_1B 
tile  diagrams  and  to  obtain  or  verify  intra-tile  offsets  Bioc[n]  and  OR 
|n]  and  0  graphically.  This  is  almost  what  has  been  done  in  Fig. 


loc  *  A12 1  and  (|A14[)  over  position-space  and  beamspace 

1  Bj)c  and  associated  index  vectors 


A1 


“Almost”  is  because  instead  of  n 
and  k  taking  values  of  0, 1,2, 3, 4,  5,  6,  in  the  figure  they  have  been  given  values  0, 1,  2, 3,  —3,  —2,  —1  for 
plotting  convenience.  This  is  acceptable  because  it  is  harmless  to  change  n  and  k  by  multiples  of  |R|.  The 
values  of  [to]  and  0  are  not  affected.  To  see  this  we  use  relationship  I  =  RR  1  =  RR,  p-|  or  its  equivalent 
RRj  =  |R|I  to  write 


nseed(n  +  m|R|)  =  nseedn  +  |R|Inseedm 
=  nseedn  +  RR*  nseedm 

=  R(m  +  Rjnseedm)  +  [to],  (A17) 


where  (|A13[)  has  been  used  in  the  last  step.  Similarly,  add  a  multiple  of  |R|  to  k  on  the  left  side  of  (|A15 1 
and  use  R*R  =  II R|  and  (IA15I)  to  obtain 


fcseed(j|R|  +  k)=  jfcseedI|R|  +  fcfcseed 
=  jfcseedRjR  +  kkseed 
=  (jfeseedRi  +  £)R  +  i 


(A  18) 


